Date   

Re: Using Python rasterio package to read gif file and world file, re-project to another picture

Sean Gillies
 

Hi Guodong,

If you save the .gif and the .gfw file to the same directory, so that they are siblings, rasterio will find the .gfw automatically when the .gif file is opened. This behavior is inherited from the GDAL library that rasterio uses and is a common GIS paradigm. Many GIS file formats are actually multi-file formats. One file, usually the image file, is the primary file and it may have auxiliary or "sidecar" files that carry additional information. A world file is one such auxiliary file.

Because rasterio's programs don't handle colormapped images well, I suggest that you use the GDAL programs instead. Like this:

gdal_translate -of GTiff -a_srs EPSG:4326 latest_radaronly.gif latest_radaronly.tif
gdalwarp -r near -t_srs EPSG:3857 latest_radaronly.tif warped.tif

Rasterio's rio-info program will show you the bounds

$ rio info warped.tif
{"bounds": [-14207635.496435506, 2470074.029222458, -7406084.60511778, 6518566.074162895], "colorinterp": ["palette"], "count": 1, "crs": "EPSG:3857", "descriptions": [null], "driver": "GTiff", "dtype": "uint8", "height": 1922, "indexes": [1], "interleave": "band", "lnglat": [-97.07967556953324, 37.395081418108624], "mask_flags": [["nodata"]], "nodata": 0.0, "res": [2106.3954448181253, 2106.3954448181253], "shape": [1922, 3229], "tiled": false, "transform": [2106.3954448181253, 0.0, -14207635.496435506, 0.0, -2106.3954448181253, 6518566.074162895, 0.0, 0.0, 1.0], "units": [null], "width": 3229}

The rio-bounds program will give you the bounds as GeoJSON

$ rio bounds warped.tif
{"bbox": [-127.62936117647057, 21.654332979481737, -66.5299899625959, 50.41561201989618], "geometry": {"coordinates": [[[-127.62936117647057, 21.654332979481737], [-66.5299899625959, 21.654332979481737], [-66.5299899625959, 50.41561201989618], [-127.62936117647057, 50.41561201989618], [-127.62936117647057, 21.654332979481737]]], "type": "Polygon"}, "properties": {"filename": "warped.tif", "id": "0", "title": "warped.tif"}, "type": "Feature"}


On Thu, Mar 5, 2020 at 6:29 AM <gd.zhu@...> wrote:

Hello,

I want to re-project a radar image (https://radar.weather.gov/Conus/RadarImg/latest_radaronly.gif) which is in NAD83/EPSG4326 to WGS84/Pseudo-Mercator/EPSG3857. There is a world file associated with this gif (https://radar.weather.gov/Conus/RadarImg/latest_radaronly.gfw). I was wondering how can I do this projection using Python rasterio package, and how to extract the coordinates bounds from the new re-projected picture?

I'm very new to GIS. I know world file contains georeferencing information, but I stucked at the first step and don't know how to read gif file along with the world file... 

Thanks a lot for your time and help! 

Best regards, 
Guodong



--
Sean Gillies


Using Python rasterio package to read gif file and world file, re-project to another picture

gd.zhu@...
 

Hello,

I want to re-project a radar image (https://radar.weather.gov/Conus/RadarImg/latest_radaronly.gif) which is in NAD83/EPSG4326 to WGS84/Pseudo-Mercator/EPSG3857. There is a world file associated with this gif (https://radar.weather.gov/Conus/RadarImg/latest_radaronly.gfw). I was wondering how can I do this projection using Python rasterio package, and how to extract the coordinates bounds from the new re-projected picture?

I'm very new to GIS. I know world file contains georeferencing information, but I stucked at the first step and don't know how to read gif file along with the world file... 

Thanks a lot for your time and help! 

Best regards, 
Guodong


Re: Creating rasterio dataset without IO

Nic Annau
 
Edited

Thanks for the suggestion.

Those functions are exactly what I would like to use, but I am having troubles with them. I think it's a rasterio issue. 

I have tried using rioxarray and XGeo, but the problem is the custom projection crs I must provide when clipping. rasterio isn't playing well with it, and I think it's worthy of opening a GitHub issue. These clipping functions seem to not like custom projections - and only work for projections with a nice EPSG reference. 

The error is: `CRSError: The PROJ4 dict could not be understood. OGR Error code 5`. The PROJ4 dict is fine so far as I can tell.

I'm purposely leaving this post sparse and will reference the GitHub issue (with more details) here. 

As for using `geometry_mask` from `rasterio.features`: It is working well! The only issues is figuring out how to define the Affine matrix without IO, and I think I'm able to do that using `rasterio.transform.from_bounds()`.

Thanks again for your time,
Nic


Re: Creating rasterio dataset without IO

Brendan Ward
 

Nic,

if your goal is to mask the data, you can replicate some of the processing steps in the `mask` chain:
https://github.com/mapbox/rasterio/blob/master/rasterio/mask.py

Specifically, you'd need to figure out the transform and output shape for your mask based on your source data and polygon, then it is just a matter of calling geometry_mask with those parameters:
https://github.com/mapbox/rasterio/blob/1.1.3/rasterio/mask.py#L108-L109

Then apply that mask to your data, which you've already prepared in a prior step.

Having the dataset as a rasterio.io.Dataset gives us access to many of the properties we need to be able to calculate those parameters.  To create a mask we don't actually need to read data from the dataset, it's only when the mask is applied that we need those data.  This use case seems outside the intent for MemoryFile.

I hope that helps.


Re: Creating rasterio dataset without IO

Alan Snow
 

Hi Nic,

I think you are looking for rioxarray. An example of what you want to do is here:
https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

Best,
Alan


Creating rasterio dataset without IO

Nic Annau
 

Hello,

The `rasterio.mask.mask` function is great for masking a raster using a polygon, however, I'm hoping to eliminate unnecessary IO and provide `rasterio.mask.mask` with a modified raster instead of loading it from file. 

Specifically, I am loading a NetCDF file as an `xarray.Dataset`, and re-gridding that file before continuing to mask it. I would love to not have to write that re-gridded `xrray` object to file before re-loading it and masking it at its new resolution.

My ultimate question is: Is there a way to convert a `numpy.ndarray`, `xarray.Dataset`, or other convenient object type into a `rasterio.io.Dataset` without having to load that data from file?

think there might be a solution somewhere using `MemoryFile`, but implementing this is not clear to me. 

My workflow is currently something like this:
```
import rasterio
import xarray as xr

ds = xr.open_dataset('path/to/file.nc')
rds = regrid_ds(ds)

# my goal is to eliminate this step
rds.to_netcdf('path/to/regridded_file.nc')

# and eliminate this step
dataset = rasterio.open('path/to/regridded_file.nc')

# and create a dataset from another object
m = mask.mask(dataset, polygon)

# continue with analysis
```

Thanks for your time!
Nic


Re: Convert gray band from dataset to RGB with rasterio for population density

Sean Gillies
 

Hi,

Sorry for the slow response.

On Fri, Feb 21, 2020 at 3:10 PM <acwangpython@...> wrote:
https://gis.stackexchange.com/questions/351614/convert-gray-band-from-dataset-to-rgb-with-rasterio-for-population-density

I'm pretty new to GIS so please correct and terminology or logical fallacies:
 
I have this dataset: https://ghsl.jrc.ec.europa.eu/download.php?ds=pop (I am using the full dataset (Global dataset) which is located in the hyperlink below the map). When I unzip it (the .tif.ovr file) and access it via rasterio, there is only one band. On QGIS GUI, I've managed to open the .tif.ovr file and change some of the colors (since I'm working on population density). Unfortunately, I have no clue on how to change this dataset to have RGB bands.
 
Right now, when I do:
 
with open(pathtodata, "r+", **profile) as src:
    src.meta
    src.dataset_mask()
 
I only get a 2D numpy array with what seems like the gray band only values (0 and 255), but I'd like to have the RGB values so I can work with the RGB values in Python (not for visualization). The meta values show that there is only one band (count) and no photometric. Doing `src.colorinterp` shows only `ColorInterp.gray: 1` which is the issue.

I downloaded a small portion of the data, GHS_POP_E2015_GLOBE_R2019A_4326_9ss_V1_0_18_4.tif, a block that covers some of the South of France and North Africa. It's the .tif file that you want to read, not the .tif.ovr file (I suspect QGIS created one for you, there are no ovr files in the downloads).

The dataset has one single float64 band. To visualize it, use a matplotlib colormap with normalization. For example:

import rasterio
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

with rasterio.open("GHS_POP_E2015_GLOBE_R2019A_4326_9ss_V1_0_18_4.tif") as dataset:
    data = src.read(1, masked=True)
    plt.imshow(data, norm=LogNorm(vmin=1e-6, vmax=data.max()))
<matplotlib.image.AxesImage object at 0x12f865710>
    plt.show()

Results in


Figure_1.png
--
Sean Gillies


Re: Convert gray band from dataset to RGB with rasterio for population density

@pyaw
 

This is a re worded question that probably sounds a bit better: https://gis.stackexchange.com/questions/351639/rasterio-change-singelband-gray-to-singleband-pseudo-color-for-multiple-colors

Quoted: This is in reference to: Convert gray band from dataset to RGB with rasterio for population density however there is a difference. I have been playing around with QGIS which I had opened the aforementioned .tif.ovr file and finally got some different hues of gray.

Instead of having multiple bands like the original question, simply having multiple colors is what I'm seeking with rasterio (or gdal. Honestly either one works). In QGIS, I set the rendering type for the .tif.ovr file to singleband pseudocolor with a magma color ramp and on continuous mode with linear interpolation. (Setting the mode to equal interval also works well).

Looks neat! So how would I do the same (probably best to use RGB instead though with black being 0 statistic minimum) with the rasterio data that I used from the aforementioned question?

Recalling, the dataset I have from src.dataset_mask() is a 2D numpy array (e.g. [[0,255],[255,255]]) in which the only values in the entire dataset is either 0 or 255. I'm seeking a way to make the values either Floats or a range between 0 and 255. Thank you! 


Re: Convert gray band from dataset to RGB with rasterio for population density

@pyaw
 

Here are some edits from stack that came about since my question may not be interpreted correctly: 

Edit: I am using the full dataset (Global dataset) which is located in the hyperlink below the map
 
Edit 2: When I mean 2D array, I meant a numpy array that looks like this: `[[0, 255, 0], [0, 0, 255]]`.
 
Additionally, this is the meta data:
`{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -200.0, 'width': 72164, 'height': 36000, 'count': 3, 'crs': None, 'transform': Affine(1.0, 0.0, 0.0, 0.0, 1.0, 0.0)}`
 
Thanks to mike for a slight correction: I am actually after just increasing the number of bands so that I can differentiate the two given values which are 0 and 255.
 
Note that when you do `x = src.dataset_mask()` to grab the numpy array, 0 and 255 are the only two values. Like any population density map, I'm after values that are between a range rather than simply having two numbers... e.g. numbers between 0-255 or float values.


Convert gray band from dataset to RGB with rasterio for population density

@pyaw
 

https://gis.stackexchange.com/questions/351614/convert-gray-band-from-dataset-to-rgb-with-rasterio-for-population-density

I'm pretty new to GIS so please correct and terminology or logical fallacies:
 
I have this dataset: https://ghsl.jrc.ec.europa.eu/download.php?ds=pop (I am using the full dataset (Global dataset) which is located in the hyperlink below the map). When I unzip it (the .tif.ovr file) and access it via rasterio, there is only one band. On QGIS GUI, I've managed to open the .tif.ovr file and change some of the colors (since I'm working on population density). Unfortunately, I have no clue on how to change this dataset to have RGB bands.
 
Right now, when I do:
 
with open(pathtodata, "r+", **profile) as src:
    src.meta
    src.dataset_mask()
 
I only get a 2D numpy array with what seems like the gray band only values (0 and 255), but I'd like to have the RGB values so I can work with the RGB values in Python (not for visualization). The meta values show that there is only one band (count) and no photometric. Doing `src.colorinterp` shows only `ColorInterp.gray: 1` which is the issue.
 
So how would I change the gray band to RGB bands to work with RGB-valued data with numpy? Thank you!


Re: Translating/shifting a raster in memory

Sean Gillies
 

Hi,

On Mon, Feb 17, 2020 at 7:20 AM <matas.lauzadis@...> wrote:
Hi all, I'm trying to apply dx, dy corrections to a DEM raster. This requires shifting the entire raster, and I was wondering if it is possible to do this in rasterio?

If you open a rasterio dataset in "r+" mode, you can modify its transform attribute. Here's a code example:

from affine import Affine
import rasterio

with rasterio.open("example.tif",  "r+") as dataset:
    original_transform = dataset.transform
    new_transform = Affine.translation(dx, dy) * original_transform
    dataset.transform = new_transform

--
Sean Gillies


Re: Using Rasterio with GDAL 2.4.x

Sean Gillies
 

Hi Christina,

On Mon, Feb 17, 2020 at 11:23 PM Ratcliff, Christina (A&F, Waite Campus) <christina.ratcliff@...> wrote:

Hi,

 

Thank you for your responses. I have previously installed rasterio via the Gohlke wheel file alongside QGIS quite successfully by using a properly configured osgeo4w shell to ensure it installs against the correct version of python and GDAL.

 

Joris, thank you for the link to QGIS on conda, but our end users aren’t programmers so installing QGIS this way may be a bit difficult for them.

 

The version of GDAL that the Rasterio wheel is built for on the Gohlke page isn’t explicitly stated, just that GDAL is required. There are wheels available for GDAL 2.2.4, 2.4.1 and 3.0.4 for python 3.7.  I have confirmed the Rasterio 1.0.24 wheel installs successfully against QGIS’s GDAL 2.4.1, 2.4.3 and 3.0.3 but Rasterio 1.0.25+ only installs against 3.0.x.

 

My method for installation is to download a wheel file from the Gohlke page and install it into the configured environment with admin privileges using

          pip install <path to rasterio wheel file>

I have attached the resulting log files created when installing Rasterio 1.0.25 & 1.1.2.

 

Just to clarify, I would like to be able to skip the building and installation of GDAL on Windows if the requirement of using GDAL >= 1.11 is already met; hence my question regarding  gdal~=3.0.1. If this is possible, where in the code is this set? Are there extra options I can use with pip to achieve this?

 

Thanks,

 

Christina


Because Gohlke's wheels are not known to pip, you must either

1) install the appropriate GDAL wheel *before* installing the rasterio wheel, or
2) use the --find-links option documented in  https://pip.pypa.io/en/stable/reference/pip_install/#cmdoption-f to point to a directory containing wheels downloaded from Gohlke's page.

If you try to install Gohlke's rasterio wheel first, pip will go looking, by default, in the Python package index for a GDAL package and will not find anything there that is compatible with Gohlke's rasterio wheel.

--
Sean Gillies


Re: Using Rasterio with GDAL 2.4.x

Ratcliff, Christina (A&F, Waite Campus)
 

Hi,

 

Thank you for your responses. I have previously installed rasterio via the Gohlke wheel file alongside QGIS quite successfully by using a properly configured osgeo4w shell to ensure it installs against the correct version of python and GDAL.

 

Joris, thank you for the link to QGIS on conda, but our end users aren’t programmers so installing QGIS this way may be a bit difficult for them.

 

The version of GDAL that the Rasterio wheel is built for on the Gohlke page isn’t explicitly stated, just that GDAL is required. There are wheels available for GDAL 2.2.4, 2.4.1 and 3.0.4 for python 3.7.  I have confirmed the Rasterio 1.0.24 wheel installs successfully against QGIS’s GDAL 2.4.1, 2.4.3 and 3.0.3 but Rasterio 1.0.25+ only installs against 3.0.x.

 

My method for installation is to download a wheel file from the Gohlke page and install it into the configured environment with admin privileges using

          pip install <path to rasterio wheel file>

I have attached the resulting log files created when installing Rasterio 1.0.25 & 1.1.2.

 

Just to clarify, I would like to be able to skip the building and installation of GDAL on Windows if the requirement of using GDAL >= 1.11 is already met; hence my question regarding  gdal~=3.0.1. If this is possible, where in the code is this set? Are there extra options I can use with pip to achieve this?

 

Thanks,

 

Christina


Translating/shifting a raster in memory

matas.lauzadis@...
 

Hi all, I'm trying to apply dx, dy corrections to a DEM raster. This requires shifting the entire raster, and I was wondering if it is possible to do this in rasterio?


Re: Using Rasterio with GDAL 2.4.x

Joris Van den Bossche
 



On Sun, 16 Feb 2020 at 01:14, Sean Gillies <sean.gillies@...> wrote:
Hi Christina,

On Thu, Feb 13, 2020 at 12:15 AM <christina.ratcliff@...> wrote:

Hi,

 

I’m trying to determine if it is possible to install Rasterio against GDAL 2.4.x on windows using the Gohlke wheels.


Yes. It is my understanding that if you download one of the GDAL wheels from https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal and then the matching rasterio wheel from https://www.lfd.uci.edu/~gohlke/pythonlibs/#rasterio and then install them using pip like `python -m pip /path/to/gdal.whl /path/to/rasterio.whl`, you can have a working rasterio installation. The rasterio wheels from that site depend on a matching gdal wheel from that same site being installed. As far as I know, Gohlke's rasterio wheels are not compatible with GDAL installed in a different way.
 

When installing rasterio versions 1.0.25 and above, it detects that the GDAL requirement is not met and attempts to build and install GDAL 3.0.4 and fails (log attached).
The environment I’m installing into is QGIS3 has GDAL 2.4.1 which shouldn’t be upgraded.


I'm pretty sure that installing one of Gohlke's rasterio wheels into your QGIS3 python will not work. I am not sure what is the best way to install rasterio into QGIS3's python on Windows. I am not even sure that it is a good idea: QGIS3 has raster capabilities and its own Python modules, it think the happy path may be to use those instead.

You can install QGIS from conda-forge nowadays, so that could be an option to install both QGIS and rasterio that should be compatible with each other (I think, didn't try on windows): https://anaconda.org/conda-forge/qgis
 
 

Is it possible to install the current release of Rasterio against GDAL 2.4.x ? The Rasterio readme suggests that this is supported.


Yes. Note that to build and install Rasterio against a GDAL version installed on your system requires headers and other development files. I am not sure if those are provided by QGIS.

On the assumption that this behavior is in error, I have been looking at source files and have a couple of queries.

 

1. In the .travis.yml are the following lines

- GDALINST=$HOME/gdalinstall
- GDALBUILD=$HOME/gdalbuild
- PROJINST=$HOME/gdalinstall
- PROJBUILD=$HOME/projbuild

 

      Should the 3rd line read $home/projinstall?


In rasterio's Travis CI configuration, PROJ and GDAL are built at different locations but then are installed into a common location which just happens to be called "gdalinstall". We could have had different install locations, but a single one is much easier to use.

 

2. I have when installing using the wheel file, that GDAL is being matched to ~=3.0.1. (see attached log)

Collecting gdal~=3.0.1

  Using cached GDAL-3.0.4.tar.gz (577 kB)


    Where is this specified in the source files? 
If Rasterio indeed supports GDAL 2.4.1 should this in fact be >= 2.x.x ?


I would need to see the command you ran to install the wheel file before I began to answer this question.

I am not familiar enough with Windows or QGIS to give really great advice. I think that QGIS users who want to use rasterio in QGIS scripts should lobby the QGIS project to start including rasterio in QGIS.
 
Many thanks in advance,


Christina


I hope I've been able to help and I hope that someone else with more Windows and QGIS knowledge will jump in to correct any misinformation that I may have given you.


Re: Using Rasterio with GDAL 2.4.x

Sean Gillies
 

Hi Christina,

On Thu, Feb 13, 2020 at 12:15 AM <christina.ratcliff@...> wrote:

Hi,

 

I’m trying to determine if it is possible to install Rasterio against GDAL 2.4.x on windows using the Gohlke wheels.


Yes. It is my understanding that if you download one of the GDAL wheels from https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal and then the matching rasterio wheel from https://www.lfd.uci.edu/~gohlke/pythonlibs/#rasterio and then install them using pip like `python -m pip /path/to/gdal.whl /path/to/rasterio.whl`, you can have a working rasterio installation. The rasterio wheels from that site depend on a matching gdal wheel from that same site being installed. As far as I know, Gohlke's rasterio wheels are not compatible with GDAL installed in a different way.
 

When installing rasterio versions 1.0.25 and above, it detects that the GDAL requirement is not met and attempts to build and install GDAL 3.0.4 and fails (log attached).
The environment I’m installing into is QGIS3 has GDAL 2.4.1 which shouldn’t be upgraded.


I'm pretty sure that installing one of Gohlke's rasterio wheels into your QGIS3 python will not work. I am not sure what is the best way to install rasterio into QGIS3's python on Windows. I am not even sure that it is a good idea: QGIS3 has raster capabilities and its own Python modules, it think the happy path may be to use those instead.
 

Is it possible to install the current release of Rasterio against GDAL 2.4.x ? The Rasterio readme suggests that this is supported.


Yes. Note that to build and install Rasterio against a GDAL version installed on your system requires headers and other development files. I am not sure if those are provided by QGIS.

On the assumption that this behavior is in error, I have been looking at source files and have a couple of queries.

 

1. In the .travis.yml are the following lines

- GDALINST=$HOME/gdalinstall
- GDALBUILD=$HOME/gdalbuild
- PROJINST=$HOME/gdalinstall
- PROJBUILD=$HOME/projbuild

 

      Should the 3rd line read $home/projinstall?


In rasterio's Travis CI configuration, PROJ and GDAL are built at different locations but then are installed into a common location which just happens to be called "gdalinstall". We could have had different install locations, but a single one is much easier to use.

 

2. I have when installing using the wheel file, that GDAL is being matched to ~=3.0.1. (see attached log)

Collecting gdal~=3.0.1

  Using cached GDAL-3.0.4.tar.gz (577 kB)


    Where is this specified in the source files? 
If Rasterio indeed supports GDAL 2.4.1 should this in fact be >= 2.x.x ?


I would need to see the command you ran to install the wheel file before I began to answer this question.

I am not familiar enough with Windows or QGIS to give really great advice. I think that QGIS users who want to use rasterio in QGIS scripts should lobby the QGIS project to start including rasterio in QGIS.
 
Many thanks in advance,


Christina


I hope I've been able to help and I hope that someone else with more Windows and QGIS knowledge will jump in to correct any misinformation that I may have given you.


Re: Trying to make custom WKT for a non-Earth projection to spatially reference images taken with microscope

Sean Gillies
 

Hi Ryan,

On Sun, Feb 9, 2020 at 1:55 PM Ryan Avery <ravery@...> wrote:

I'm working with non spatially-referenced images produced from a laser ablation system. Along with the image, the software produces an "align" file of the same name - an xml file that contains the XY center point of each image (coordinates are in microns) and the XY size (in microns) along with any rotation.

I would like to display the images (we typically generate about 1000 images / day) in their correct relative positions and sizes.

The laser has a stage with a cartesian positive XY coordinate system (in microns), and the origin in the SW corner.

So far I have gotten to the step of generating an affine transform for each scanned image based on it's size and center coordinate in microns and rotation metadata. What I don't have is a definition of a custom CRS in WKT. I don't have any experience with the WKT format but it seems like the way to represent non-Earth CRSs. Once I have this I think I reproject each image and save as a geotiff for later viewing in QGIS.

Are there any resources for defining a custom CRS in WKT or has anyone done something like this before? I know the units are in microns, the center point of my CRS, and the height and width of the coordinate reference system (and that it is planar). Any help is much appreciated, my progress so far is located here: https://github.com/rbavery/custom_spatial_reference/blob/master/spatial_ref_imgs.ipynb

Best,
Ryan


This sounds like a job for WKT2 EngineeringCRS.


I do not know how to use it with rasterio, but I would love to learn how.

--
Sean Gillies


Using Rasterio with GDAL 2.4.x

Ratcliff, Christina (A&F, Waite Campus)
 

Hi,

 

I’m trying to determine if it is possible to install Rasterio against GDAL 2.4.x on windows using the Gohlke wheels.

 

When installing rasterio versions 1.0.25 and above, it detects that the GDAL requirement is not met and attempts to build and install GDAL 3.0.4 and fails (log attached).
The environment I’m installing into is QGIS3 has GDAL 2.4.1 which shouldn’t be upgraded.

 

Is it possible to install the current release of Rasterio against GDAL 2.4.x ? The Rasterio readme suggests that this is supported.

 

On the assumption that this behavior is in error, I have been looking at source files and have a couple of queries.

 

1. In the .travis.yml are the following lines

- GDALINST=$HOME/gdalinstall
- GDALBUILD=$HOME/gdalbuild
- PROJINST=$HOME/gdalinstall
- PROJBUILD=$HOME/projbuild

 

      Should the 3rd line read $home/projinstall?

 

2. I have when installing using the wheel file, that GDAL is being matched to ~=3.0.1. (see attached log)

Collecting gdal~=3.0.1

  Using cached GDAL-3.0.4.tar.gz (577 kB)


    Where is this specified in the source files? 
If Rasterio indeed supports GDAL 2.4.1 should this in fact be >= 2.x.x ?


Many thanks in advance,

Christina


Trying to make custom WKT for a non-Earth projection to spatially reference images taken with microscope

Ryan Avery
 

I'm working with non spatially-referenced images produced from a laser ablation system. Along with the image, the software produces an "align" file of the same name - an xml file that contains the XY center point of each image (coordinates are in microns) and the XY size (in microns) along with any rotation.

I would like to display the images (we typically generate about 1000 images / day) in their correct relative positions and sizes.

The laser has a stage with a cartesian positive XY coordinate system (in microns), and the origin in the SW corner.

So far I have gotten to the step of generating an affine transform for each scanned image based on it's size and center coordinate in microns and rotation metadata. What I don't have is a definition of a custom CRS in WKT. I don't have any experience with the WKT format but it seems like the way to represent non-Earth CRSs. Once I have this I think I reproject each image and save as a geotiff for later viewing in QGIS.

Are there any resources for defining a custom CRS in WKT or has anyone done something like this before? I know the units are in microns, the center point of my CRS, and the height and width of the coordinate reference system (and that it is planar). Any help is much appreciated, my progress so far is located here: https://github.com/rbavery/custom_spatial_reference/blob/master/spatial_ref_imgs.ipynb

Best,
Ryan


Re: Unexpected result using rasterio windowed writing - Access window out of range in RasterIO()

umbertofilippo@...
 

Replying to myself to show my last update on the situation...

I have made some interesting improvements:

I am now able to write the full extent of raster `A` correctly over raster `B`, having the same pixel values in the same locations for both rasters.

The only thing left before having a complete solution is that the output raster has some (apparently regularly spaced) stripes, sign that there might be some adjustments to do in the window coordinates translation.

Here is a picture of the updated situation:



zoomed to see the stripes:



Finally, the updated code:

    with rasterio.open(inputfile, 'r') as src:
        with rasterio.open(
            os.path.join(dest_dir, 'mosaic2.tif'),
            'w',
            driver='GTiff',
            compress="DEFLATE",
            height=56160,
            width=85196,
            count=1,
            dtype=np.float64,
            crs='+proj=latlong',
            transform=out_transform,
        ) as mosaic_raster:
            for ji, window in src.block_windows(1):
                r = src.read(1, window=window)
                left, top = src.xy(
                    row=window.row_off,
                    col=window.col_off)
                right, bottom = src.xy(
                    row=window.row_off + window.height,
                    col=window.col_off + window.width)
                new_window = from_bounds(
                    left=left,
                    bottom=bottom,
                    right=right,
                    top=top,
                    transform=mosaic_raster.transform,
                    height=window.height,
                    width=window.width)
                mosaic_raster.write(r, 1, window=new_window)