Date   

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


Re: Interconverting between pyproj and rasterio CRS objects

Alan Snow
 

Here are instructions for rasterio to geopandas: https://gis.stackexchange.com/a/350231/144357

For geopandas to rasterio, you will need something like:

 
if LooseVersion(rasterio.__gdal_version__) > LooseVersion("3.0.0"):
 
crs = rasterio.crs.CRS.from_wkt(crs.to_wkt())
else: 
crs = rasterio.crs.CRS.from_wkt(crs.to_wkt("WKT1_GDAL"))


Apologies for the formatting ....


Interconverting between pyproj and rasterio CRS objects

Nick Weir
 

Hi there!

Now that geopandas is using pyproj CRS objects as the default for tracking CRSs in GeoDataFrames, I'm finding myself needing to figure out a way to frequently check identity and interconvert back and forth between pyproj.crs.CRS objects and rasterio.crs.CRS objects. However, I'm having a hard time figuring out the "ideal" way to do this. It's great that the following works:
import rasterio
import pyproj
pp_crs = pyproj.crs.CRS.from_epsg(4326)
rio_crs = rasterio.crs.CRS.from_epsg(4326)
assert pp_crs == rio_crs

I'm trying to figure out the best way to pass around CRSs within a python library (ideally in a single object type) and encountering issues. I don't want to wed myself to using EPSG codes to interconvert between objects, because I want to be able to accommodate CRSs/SRSs without codes. However, I'm finding problems using something else (like a WKT string) to convert between rasterio.crs.CRSs and pyproj.crs.CRSs. For example:
import rasterio
import pyproj
pp_crs = pyproj.crs.CRS.from_epsg(4326)
rio_crs = rasterio.crs.CRS.from_epsg(4326)
pp_to_rio_crs = rasterio.crs.CRS.from_wkt(pp_crs.to_wkt())
assert rio_crs == pp_to_rio_crs
Fails. Any thoughts? I could definitely be missing something here. Thanks! -N


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.