Re: New rasterio versions don't handle blank SRS gracefully

Sean Gillies

Hi Nicholas,

Thanks for the explanation. We have at least one more related regression in 1.0.15 and I intend to fix both of these in 1.0.16.

On Fri, Feb 1, 2019 at 10:36 AM <nicholas.maxwell@...> wrote:
Short version:

The following causes errors, in rasterio 1.0.15,  but is fine in rasterio 1.0.9. 

from import CRS
The following are the only tests in rasterio to check blank CRSs.

def test_null_crs_equality():
assert CRS() == CRS()
a = CRS()
assert a == a
assert not a != a

def test_null_and_valid_crs_equality():
assert (CRS() == CRS(init='epsg:4326')) is False

I feel like print(CRS()) should not cause errors - CRS().__str__ should just return an empty string, for example.

Longer version:

We work with pre-orthorectified sat imagery, which have blank geotransforms (but may have GCPs set), eg,

gdalinfo /data/16APR27183709-P1BS-502376654080_01_P009.tif 
Driver: GTiff/GeoTIFF
Files: /data/16APR27183709-P1BS-502376654080_01_P009.tif
Size is 3201, 3281
Coordinate System is `'
GCP Projection = 
        SPHEROID["WGS 84",6378137,298.257223563,

This causes problems on newer rasterio versions (presumably due to the switch to WKT representation for CRS). Example:

import rasterio
with'/data/16APR27183709-P1BS-502376654080_01_P009.tif') as src:

On rasterio 1.0.9, this is fine, evaluates to None:

python ~/
/home/nick/venvs/ras9/lib/python3.6/site-packages/rasterio-1.0.9-py3.6-linux-x86_64.egg/rasterio/ NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: None

The NotGeoreferencedWarning is fine, we explicitly ignore it. The point is that is None.
On 1.0.15, however,

python ~/
/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/ NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.
  s = DatasetReader(path, driver=driver, **kwargs)
CRS: Traceback (most recent call last):
  File "rasterio/_crs.pyx", line 274, in rasterio._crs._CRS.to_dict
  File "rasterio/_err.pyx", line 194, in rasterio._err.exc_wrap_ogrerr
rasterio._err.CPLE_BaseError: OGR Error code 7

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/nick/", line 3, in <module>
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/", line 245, in to_string
    return self.to_wkt() or self.to_proj4()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/", line 106, in to_proj4
    return ' '.join(['+{}={}'.format(key, val) for key, val in])
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/", line 170, in data
    self._data = self.to_dict()
  File "/home/nick/venvs/ras15/lib/python3.6/site-packages/rasterio/", line 164, in to_dict
    return self._crs.to_dict()
  File "rasterio/_crs.pyx", line 277, in rasterio._crs._CRS.to_dict
rasterio.errors.CRSError: The WKT could not be parsed. OGR Error code 7

Interestingly, the GDAL Python API handles this in a different way:

from osgeo import gdal
ds = gdal.Open('/data/16APR27183709-P1BS-502376654080_01_P009.tif')
prj = ds.GetProjection()

which prints

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]
I guess GDAL assumes the GCP projection, in case of blank geotransform.

In any case, I think the is None behavior is good, since we can just look at src.gcps for the gcps.

I will add a workaround on our side, since
works, (that is, CRS().is_valid evaluates to False, without throwing exceptions.

Sean Gillies

Join to automatically receive all group messages.