I'm trying to write the output of some raster processing using rasterio 0.36.0 (conda package, Python 2.7, Windows 32 bit). The output raster should have the same metadata as the input, including the coordinate reference system (crs). However when I set the output metadata using **kwargs, the crs is missing from the output.

See the simplified example below (without any further processing) importing then exporting a subset of a Landsat8 blue band:

import rasterio

raster_in = "Blueband_crestmap_sample.tif"
raster_out = "Blue_sample_copy.tif"

with rasterio.open(raster_in) as inrast:
    blue = inrast.read(1)
    kwargs = inrast.profile
    print inrast.crs
    print kwargs

with rasterio.open(raster_out, 'w', **kwargs) as outrast:
    outrast.write_band(1, blue)
    print outrast.crs
    print outrast.profile

Input dataset:

CRS({'init': u'epsg:32653'})

{'count': 1, 'crs': CRS({'init': u'epsg:32653'}), 'interleave': 'band', 'dtype': 'int16', 'affine': Affine(30.0, 0.0, 648855.0, 0.0, -30.0, -2714235.0), 'driver': u'GTiff', 'transform': (648855.0, 30.0, 0.0, -2714235.0, 0.0, -30.0), 'height': 440, 'width': 503, 'tiled': False, 'nodata': -9999.0}

Output dataset:


{'count': 1, 'crs': CRS({}), u'interleave': 'band', 'dtype': 'int16', 'affine': Affine(30.0, 0.0, 648855.0, 0.0, -30.0, -2714235.0), 'driver': u'GTiff', 'transform': (648855.0, 30.0, 0.0, -2714235.0, 0.0, -30.0), 'height': 440, 'width': 503, u'tiled': False, 'nodata': -9999.0}

Aside from the crs, the writing is successful. For example, I can very easily set the crs from the catalog in ArcGIS, but this is not a good solution for multiple datasets.

This is not specific to Landsat. I have the same issue when running a filter to smooth a DEM.

Is there anything I'm missing? The code seems straight forward enough and matches a similar workflow from the rasterio cookbook: https://mapbox.s3.amazonaws.com/playground/perrygeo/rasterio-docs/cookbook.html#using-a-scipy-filter-to-smooth-a-raster

  • The CRS needs to be set separately, use GetProjection / SetProjection like in gis.stackexchange.com/questions/60371/… – Michael Stimson Apr 3 '18 at 5:25
  • Thanks @MichaelStimson. GDAL can be used to set the projection as you suggest. However, rasterio should be able to achieve simple tasks like this without calling having to directly call GDAL in my script. I managed to solve the issue. See answer below. – user96805 Apr 4 '18 at 6:24

I found a related discussion here.

I was able to solve the issue by setting my GDAL_DATA environment variable, so that rasterio (gdal) could find the gcs.csv file containing epsg info. See discussion here.

It seems to be an issue that arises occasionally when installing gdal with conda on Windows (reported here).

  • Are you installing rasterio and running your code in your root conda env? rasterio pulls in gdal which installs an init script that sets the required gdal environment variables when a conda env is activated but the root conda doesn't run these init scripts. – user2856 Apr 6 '18 at 12:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.