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:
CRS({})
{'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
