Date   

Re: Any way to ensure cell resolutions stay square?

Luke
 

Try going the other way - multiply the resolution by a scale factor (i.e. 3 in your case) and divide the dimensions - see example gist

On Fri, 13 Mar 2020 at 14:35, himat15 via Groups.Io <himat15=yahoo.com@groups.io> wrote:
 
I am resampling a raster to become a different resolution, but it seems like the resolutions are not entirely accurate.

Here is the code I have (the code is very standard, taken from the documentation).

orig_res = ds.res[0] # Calculate scale factor for this specific input layer upscale_factor = orig_res / 15.0 print("upscale factor: ", upscale_factor) # Scale while reading in array scaled_arr = ds.read( out_shape=( ds.count, int(ds.height * upscale_factor), int(ds.width * upscale_factor) ), resampling=Resampling.bilinear ) # Also have to update the transform to be scaled scaled_transform = ds.transform * ds.transform.scale( (ds.height / scaled_arr.shape[1]), (ds.width / scaled_arr.shape[2]) ) scaled_profile = ds.profile scaled_profile.update(transform=scaled_transform, height=scaled_arr.shape[1], width=scaled_arr.shape[2]) with rio.open(f_out_path, 'w', **scaled_profile) as scaled_ds: scaled_ds.write(scaled_arr) print(f"Orig arr shape: ({ds.count}, {ds.shape[0]}, {ds.shape[1]})") print(f"Scaled ds shape: ({scaled_ds.count}, {scaled_ds.shape[0]}, {scaled_ds.shape[1]})") print("Orig ds res: ", ds.res) print("Scaled ds res: ", scaled_ds.res)

Code output:
upscale factor:  0.3333333333333333
Orig arr shape: (1, 9240, 2671)
Scaled ds shape: (1, 3080, 890)
Orig ds res:  (5.0, 5.0)
Scaled ds res:  (15.0, 15.00561797752809)

But in the final scaled ds resolution, why is the vertical scale direction slightly different? (15.0, 15.00561797752809)?

I'm not sure if it's a bad thing or not. I just want to make sure it's okay, or if there's a way to correct it, that would be good too.


Any way to ensure cell resolutions stay square?

himat15@...
 

 
I am resampling a raster to become a different resolution, but it seems like the resolutions are not entirely accurate.

Here is the code I have (the code is very standard, taken from the documentation).

orig_res = ds.res[0] # Calculate scale factor for this specific input layer upscale_factor = orig_res / 15.0 print("upscale factor: ", upscale_factor) # Scale while reading in array scaled_arr = ds.read( out_shape=( ds.count, int(ds.height * upscale_factor), int(ds.width * upscale_factor) ), resampling=Resampling.bilinear ) # Also have to update the transform to be scaled scaled_transform = ds.transform * ds.transform.scale( (ds.height / scaled_arr.shape[1]), (ds.width / scaled_arr.shape[2]) ) scaled_profile = ds.profile scaled_profile.update(transform=scaled_transform, height=scaled_arr.shape[1], width=scaled_arr.shape[2]) with rio.open(f_out_path, 'w', **scaled_profile) as scaled_ds: scaled_ds.write(scaled_arr) print(f"Orig arr shape: ({ds.count}, {ds.shape[0]}, {ds.shape[1]})") print(f"Scaled ds shape: ({scaled_ds.count}, {scaled_ds.shape[0]}, {scaled_ds.shape[1]})") print("Orig ds res: ", ds.res) print("Scaled ds res: ", scaled_ds.res)

Code output:
upscale factor:  0.3333333333333333
Orig arr shape: (1, 9240, 2671)
Scaled ds shape: (1, 3080, 890)
Orig ds res:  (5.0, 5.0)
Scaled ds res:  (15.0, 15.00561797752809)

But in the final scaled ds resolution, why is the vertical scale direction slightly different? (15.0, 15.00561797752809)?

I'm not sure if it's a bad thing or not. I just want to make sure it's okay, or if there's a way to correct it, that would be good too.


Re: Decimated and windowed read in rasterio

Sean Gillies
 

Hi,

On Thu, Mar 12, 2020 at 9:30 AM <umbertofilippo@...> wrote:

(Posted also on gis.stackexchange.com: https://gis.stackexchange.com/questions/353794/decimated-and-windowed-read-in-rasterio)

In my workflow using rasterio, I'd like to read an overview from a raster and get only a portion of it through a window at the same time.

Is this possible?

I have a pretty complex script so far and I'm trying to check whether the wrong output is due to this not being possible.

Basically, I am asking if it is possible to do something like this:

data = src.read(out_shape=(1,
                           math.ceil(src.height / 64),
                           math.ceil(src.width / 64)),
                window=window)

where data is a decimated read of src to get source raster overview for factor 64, and window has the width and height whose dimensions are relative to the source raster (not the overview).


the out_shape and window parameters of a Rasterio dataset's read method are very similar to the nBufXSize, nBufYSize, nXOff, nYOff, nXSize, and nYSize parameters of GDAL's RasterIO function. The window parameters refer to a region of the source raster at its full resolution, overviews are disregarded. The out_shape paramter determines the size of the output buffer and GDAL *may* read from overviews when they are available to reduce the cost of resampling.

Rasterio's read(), like GDAL's RasterIO, abstracts the details of full resolution vs overviews. Depending on the math, you may or may not get exactly the overview level you are looking for. If you absolutely need a particular overview level, I suggest passing an overview_level keyword argument to rasterio.open().

--
Sean Gillies


Decimated and windowed read in rasterio

umbertofilippo@...
 

(Posted also on gis.stackexchange.com: https://gis.stackexchange.com/questions/353794/decimated-and-windowed-read-in-rasterio)

In my workflow using rasterio, I'd like to read an overview from a raster and get only a portion of it through a window at the same time.

Is this possible?

I have a pretty complex script so far and I'm trying to check whether the wrong output is due to this not being possible.

Basically, I am asking if it is possible to do something like this:

data = src.read(out_shape=(1,
                           math.ceil(src.height / 64),
                           math.ceil(src.width / 64)),
                window=window)

where data is a decimated read of src to get source raster overview for factor 64, and window has the width and height whose dimensions are relative to the source raster (not the overview).


Rasterio Python Library Merging R, G, B Bands to Create a GeoReferenced GeoTiff

nathan.raley@...
 

I am attempting to use rasterio to merge 3 bands from Sentinel-2 data, sourced as .JP2 files, into a georeferenced GeoTiff.  However, I am running into 2 problems doing this.

First of all, here is my code for the merge:
    #Setup the bands
    b4 = rio.open(path + Bands.B04.value[0])
    b3 = rio.open(path + Bands.B03.value[0])
    b2 = rio.open(path + Bands.B02.value[0])
   
    #Setup the path for the constructed RGB image
    tif = path + "RGB.tif"
   
    #Create the RGB image
    with rio.open(tif, 'w', driver='Gtiff', width=b4.width, height=b4.height, count=3, crs=b4.crs, tranform=b4.transform, dtype=b4.dtypes[0], photometric="RGB", tile=True, compress='lzw', bigtiff=True) as rgb:
        rgb.write(b2.read(1), 1)
        rgb.write(b3.read(1), 2)
        rgb.write(b4.read(1), 3)
        rgb.close()

Now, how do I copy the source projection information and bounds in the new GeoTiff I am creating?  Second, why is it showing up as all black when viewing the image in most common image browsers, but just fine in mapping software such as ArcMap or ArcGIS Pro?

Any ideas?


Re: Transverse Mercator Complex projection

Luke
 

No. Exactly as I show in my code example. No need to write it out again, as you are just updating the georeferencing with the dataset opened in update mode ("r+").  The dataset is saved automatically when you exit the `with` statement context.

Luke


Re: Do I need to warp jp2 satellite images for rasterio to manipulate them?

himat15@...
 

What is the best way to do this warping?

Should I not actually warp and just set the `.crs` property of the jp2 raster and then re-save it in the same file?
Or should I warp it the crs, and then re-save as a .jp2 or should I save as a .tif?

Should I use virtual warping like here https://rasterio.readthedocs.io/en/latest/topics/virtual-warping.html
or should I use reprojection like here https://rasterio.readthedocs.io/en/latest/topics/reproject.html 

I just want to make sure I don't mess anything up, so thank you!


Re: Transverse Mercator Complex projection

himat15@...
 

@ Luke Pinner,

How should I then save the file with the updated CRS?
Like this? Or is there a better way?

with rasterio.open(in_tif, 'r+') as src:
src.crs = CRS.from_epsg(32629)
with raster.open("out.tif", 'w', **src.profile) as dst:
dst.write(src.read(1))


Re: Transverse Mercator Complex projection

Luke
 

To override the projection with rasterio:

import rasterio
from rasterio.crs import CRS

in_tif = 'path/to.tif'
with rasterio.open(in_tif, 'r+') as src:
src.crs = CRS.from_epsg(32629)


Re: Transverse Mercator Complex projection

Even Rouault
 

On lundi 9 mars 2020 14:42:32 CET himat15 via Groups.Io wrote:
Hi, I have a geotiff with a "Transverse_Mercator_Complex" CRS.

Here is the output of the `.crs` property after reading in the tif with
rio.read():
PROJCS["WGS_1984_Complex_UTM_Zone_29N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_
1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0
],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator_Compl
ex"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],P
ARAMETER["Central_Meridian",-9.0],PARAMETER["Scale_Factor",0.9996],PARAMET
ER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
The rest of my files have the crs EPSG:32629, so I tried to reproject this
tif into EPSG:32629 following the standard given example in the docs here (
https://rasterio.readthedocs.io/en/latest/topics/reproject.html ).

But, I got this error:
No translation for Transverse_Mercator_Complex to PROJ.4 format is known

I found that this projection is the same (?) as the regular Transverse
Mercator from
https://gis.stackexchange.com/questions/226679/complex-utm-projection. (
https://gis.stackexchange.com/questions/226679/complex-utm-projection ) It
seems to just be an ESRI specific one. It's also listed in this repo
https://github.com/Esri/projection-engine-db-doc/blob/master/text/pe_list_p
rojection.txt. (
https://github.com/Esri/projection-engine-db-doc/blob/master/text/pe_list_p
rojection.txt )

What should I do in this case to transform my file in this case
to EPSG:32629?
My understanding is that Transverse_Mercator_Complex corresponds to PROJ
"etmerc" method, which is the default since PROJ 4.9.3 for UTM. PROJ just
lacks a mapping for the Transverse_Mercator_Complex term.

So you could just override the projection information with for example:
gdal_translate in.tif out.tif -a_srs EPSG:32629
(or gdal_edit.py in.tif -a_srs EPSG:32629)
There's no need to do a reprojection.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com


Transverse Mercator Complex projection

himat15@...
 

Hi, I have a geotiff with a "Transverse_Mercator_Complex" CRS.

Here is the output of the `.crs` property after reading in the tif with rio.read():
PROJCS["WGS_1984_Complex_UTM_Zone_29N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator_Complex"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]
The rest of my files have the crs EPSG:32629, so I tried to reproject this tif into EPSG:32629 following the standard given example in the docs here.

But, I got this error:
No translation for Transverse_Mercator_Complex to PROJ.4 format is known

I found that this projection is the same (?) as the regular Transverse Mercator from
https://gis.stackexchange.com/questions/226679/complex-utm-projection. It seems to just be an ESRI specific one. It's also listed in this repo https://github.com/Esri/projection-engine-db-doc/blob/master/text/pe_list_projection.txt.

What should I do in this case to transform my file in this case to EPSG:32629?


Re: Do I need to warp jp2 satellite images for rasterio to manipulate them?

Sean Gillies
 

Hi,

On Mon, Mar 9, 2020 at 2:05 PM himat15 via Groups.Io <himat15=yahoo.com@groups.io> wrote:
I have some .jp2 files that don't have a `.crs` property after I used `rio.open` on those files. 
I read https://gis.stackexchange.com/a/233417/143989 and after using the `.gcps` property, I was able to get the crs.

But, the second answer also suggested that I could use `gdalwarp` to embed the projection info in the raster directly.
My question is, do I need to do the `gdalwarp`? I don't know if not having the crs in the normal property will affect rasterio's ability to do other operations on the jp2 file raster such as transformations to change the cell size. Or will all rasterio operations work as normal with the jp2 file as it is?

I recommend warping. GCPs are not seamlessly supported throughout rasterio 1.1.3. For example, a dataset's window_bounds and window_transform methods are ignorant of the dataset's GCPs and their CRS.

--
Sean Gillies


Do I need to warp jp2 satellite images for rasterio to manipulate them?

himat15@...
 

I have some .jp2 files that don't have a `.crs` property after I used `rio.open` on those files. 
I read https://gis.stackexchange.com/a/233417/143989 and after using the `.gcps` property, I was able to get the crs.

But, the second answer also suggested that I could use `gdalwarp` to embed the projection info in the raster directly.
My question is, do I need to do the `gdalwarp`? I don't know if not having the crs in the normal property will affect rasterio's ability to do other operations on the jp2 file raster such as transformations to change the cell size. Or will all rasterio operations work as normal with the jp2 file as it is?



Re: Same file doesn't work on non-s3 host

Sean Gillies
 

Hi Olivier,

On Thu, Mar 5, 2020 at 10:07 AM Olivier Dalang <olivier.dalang@...> wrote:
Dear list,

I'm just starting with rasterio (through the marblecutter project), and am stumbling on a strange issue. Everything works fine when I try to load tifs from s3 (I'm testing out with openaerialmap's data). But if I serve one of their files from another source, it doesn't work, with the following error :

rasterio.errors.RasterioIOError: '/vsicurl/https://[...].tif' not recognized as a supported file format.
I looked carefully at the data and headers returned by both s3 and my server, and I don't see any significant difference, neither in the data or the headers...

Is there some s3 specific logic in rasterio that prevents this from working on a regular http server ?

Or how can I get some more info about what is failing ? I'm a bit lost about how to debug this...

Thanks !!

Olivier

The exception message doesn't give you much of a lead, it's true. To debug network issues, the places to look are your server logs and GDAL's curl client logs. You can enable the latter by running your program with CPL_CURL_VERBOSE=YES set in the environment. Summaries of the HTTP requests and responses will dumped into your terminal, as if you ran "curl -v".

--
Sean Gillies


Re: Interconverting between pyproj and rasterio CRS objects

Nick Weir
 

Got it. Forgot that there wasn’t a class there previously...thanks for the context!
--
Sent from my mobile.


Re: Interconverting between pyproj and rasterio CRS objects

Sean Gillies
 

Hi Nick,

Rasterio has its own CRS class because pyproj 2 didn't exist at the time. 

On Thu, Mar 5, 2020 at 9:47 AM Nick Weir <nicholas.r.weir@...> wrote:
Yeah, this of course works. I'm more thinking from a philosophical standpoint of "which object works better to pass around?"


It seems like the pyproj objects work a bit better. And that's what's natively used in geopandas and a bunch of other places. I guess the question becomes "what's the benefit of rasterio having its own CRS object type"...I guess the main argument against is that adding pyproj is adding an additional dependency?

See this gist for a bit of exploration of this.



--
Sean Gillies


Re: Interconverting between pyproj and rasterio CRS objects

Alan Snow
 

Add some history for reference: https://rasterio.groups.io/g/dev/topic/switch_from_proj_4_format_to/28991934


Same file doesn't work on non-s3 host

olivier.dalang@...
 

Dear list,

I'm just starting with rasterio (through the marblecutter project), and am stumbling on a strange issue. Everything works fine when I try to load tifs from s3 (I'm testing out with openaerialmap's data). But if I serve one of their files from another source, it doesn't work, with the following error :

rasterio.errors.RasterioIOError: '/vsicurl/https://[...].tif' not recognized as a supported file format.
I looked carefully at the data and headers returned by both s3 and my server, and I don't see any significant difference, neither in the data or the headers...

Is there some s3 specific logic in rasterio that prevents this from working on a regular http server ?

Or how can I get some more info about what is failing ? I'm a bit lost about how to debug this...

Thanks !!

Olivier


Re: Interconverting between pyproj and rasterio CRS objects

Nick Weir
 

Yeah, this of course works. I'm more thinking from a philosophical standpoint of "which object works better to pass around?"


It seems like the pyproj objects work a bit better. And that's what's natively used in geopandas and a bunch of other places. I guess the question becomes "what's the benefit of rasterio having its own CRS object type"...I guess the main argument against is that adding pyproj is adding an additional dependency?

See this gist for a bit of exploration of this.


Re: Interconverting between pyproj and rasterio CRS objects

vincent.sarago@...
 

FYI: https://github.com/corteva/rioxarray/issues/92 and https://github.com/corteva/rioxarray/issues/92

241 - 260 of 698