1. Rasterio | Main
  2. Topics

Warp doesn't work with a general oblique transform?

erik.koene@...
 

The following Python code...
import rasterio
from rasterio.crs import CRS
from rasterio import shutil as rio_shutil
from rasterio.vrt import WarpedVRT
import matplotlib.pyplot as plt

path = 'CORINE/U2018_CLC2018_V2020_20u1.tif'

# Destination CRS is 4326
# dst_crs = CRS.from_epsg(4326) # <--- This one will work...
dst_crs = CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137') # <--- This one doesn't work

# Output image dimensions
dst_height = dst_width = 1000

# Output image transform
left, bottom, right, top = 10.897971,49.059812, 10.907593,49.064142 # Coordinates in target CRS
xres = (right - left) / dst_width
yres = (top - bottom) / dst_height
dst_transform = rasterio.Affine(xres, 0.0, left,
                              0.0, -yres, top)

vrt_options = {
    'crs': dst_crs,
    'transform': dst_transform,
    'height': dst_height,
    'width': dst_width,
}

with rasterio.open(path) as src:
    with WarpedVRT(src, **vrt_options) as vrt:
        # Read all data into memory and show plot
        data = vrt.read()
        plt.imshow(data.squeeze())
fails to run, giving the error message
The PROJ4 dict could not be understood. OGR Error code 5
which doesn't appear if I use a more 'standard' destination CRS, e.g., CRS.from_epsg(4326) . I'll welcome any suggestions! Thanks.
Sean Gillies
 

Hi,

With the latest release of GDAL and PROJ 8 I see

>>> CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137')
CRS.from_wkt('GEOGCRS["unnamed",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],DERIVINGCONVERSION["unknown",METHOD["PROJ ob_tran o_proj=longlat"],PARAMETER["lon_0",-170,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lon_p",-180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lat_p",137,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]]')

What versions of rasterio, GDAL, and PROJ do you have? And where did you install them from?


On Wed, Oct 13, 2021 at 10:04 AM <erik.koene@...> wrote:
The following Python code...
import rasterio
from rasterio.crs import CRS
from rasterio import shutil as rio_shutil
from rasterio.vrt import WarpedVRT
import matplotlib.pyplot as plt

path = 'CORINE/U2018_CLC2018_V2020_20u1.tif'

# Destination CRS is 4326
# dst_crs = CRS.from_epsg(4326) # <--- This one will work...
dst_crs = CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137') # <--- This one doesn't work

# Output image dimensions
dst_height = dst_width = 1000

# Output image transform
left, bottom, right, top = 10.897971,49.059812, 10.907593,49.064142 # Coordinates in target CRS
xres = (right - left) / dst_width
yres = (top - bottom) / dst_height
dst_transform = rasterio.Affine(xres, 0.0, left,
                              0.0, -yres, top)

vrt_options = {
    'crs': dst_crs,
    'transform': dst_transform,
    'height': dst_height,
    'width': dst_width,
}

with rasterio.open(path) as src:
    with WarpedVRT(src, **vrt_options) as vrt:
        # Read all data into memory and show plot
        data = vrt.read()
        plt.imshow(data.squeeze())
fails to run, giving the error message
The PROJ4 dict could not be understood. OGR Error code 5
which doesn't appear if I use a more 'standard' destination CRS, e.g., CRS.from_epsg(4326) . I'll welcome any suggestions! Thanks.



--
Sean Gillies
erik.koene@...
 

Hi Sean, thanks for your response. I'm running Python 3.7, rasterio version 1.0.14; gdalinfo reports version 2.4, and proj reports version 5.2.0. I'd like to clarify with you that this is a result from rasterio.crs.CRS.from_string, and not from another library? Because when I use pyproj, the CRS is in principle created correctly,

So, question 1: how do you generate the CRS, with pyproj or with rasterio? Below, you can see that the pyproj call returns no problem yet, while the rasterio version already does.
import pyproj
dst_crs = pyproj.CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137')
print(dst_crs) # >> +proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137 +type=crs
import rasterio
dst_crs = rasterio.crs.CRS.from_string('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137')

Question 2, when I use rasterio.crs.CRS.from_wkt() and use the string you post, then the warp still returns an error:
---------------------------------------------------------------------------
CRSError                                  Traceback (most recent call last)
<ipython-input-6-8658800cad11> in <module>
     29 
     30 with rasterio.open(path) as src:
---> 31     with WarpedVRT(src, **vrt_options) as vrt:
     32         # Read all data into memory and show plot
     33         data = vrt.read()

rasterio/_warp.pyx in rasterio._warp.WarpedVRTReaderBase.__init__()

rasterio/_base.pyx in rasterio._base._osr_from_crs()

CRSError: Invalid CRS: CRS.from_wkt('GEOGCRS["unnamed",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],DERIVINGCONVERSION["unknown",METHOD["PROJ ob_tran o_proj=longlat"],PARAMETER["lon_0",-170,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lon_p",-180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],PARAMETER["o_lat_p",137,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CS["ellipso
It complains about an "invalid" CRS at the stage of the virtual warp, even though the WKT string is accepted.

The oblique transform works fine to transform individual coordinates. I.e., running the following snippet
from pyproj import Transformer
rotated_to_CRS = Transformer.from_crs('+proj=ob_tran +o_proj=longlat +lon_0=-170 +o_lon_p=-180 +o_lat_p=137',
        'EPSG:4326',
        always_xy=True)
coordinates=[(0,0), (1,1)]
[pt for pt in rotated_to_CRS.itransform(coordinates)] # >> [(9.99999999999999, 47.0), (11.49406628349866, 47.990464430559726)]
correctly transforms from the oblique rotated coordinate system into the WGS84 system. So I don't understand what is strictly wrong with the CRS.
Denis Rykov
 

I would recommend you to upgrade rasterio. The version you are using is quite old.
erik.koene@...
 

Hi Denis, I'm ashamed and happy to say that that indeed did the trick, the warp works perfectly now. Thanks!

1 - 5 of 5