gyanz/pydsstools

readme has incorrect example for overriding CRS when calling `save_tiff`

Closed this issue · 3 comments

mxkpp commented

dataset.raster.save_tiff(r'grid_dataset.tif', crs=2868)

Should be like

dataset.raster._obj._crs = pyproj.CRS.from_epsg(2868)
dataset.raster.save_tiff(r'grid_dataset.tif')
mxkpp commented

I just noticed there's an open PR to update this part of the readme.

Please note that for DSS grids whose original CRS fails in rasterio to cast to wkt, the whole situation fails before it gets a chance to parse the custom profile dictionary provided to the save_tiff method. See screenshot of traceback below.

image

In this case, the approach I suggested above, modifying dataset.raster._obj._crs before calling save_tiff, is the way to go.

Agreed. My correction would not fix if the DSS grid CRS code fails to parse to rasterio CRS wkt. My correction was only to fix the TypeError since the second argument for save_tiff was a rasterio profile rather than the crs code itself.

I was utilizing a more elaborate way of fixing this issue by not even calling the save_tiff method.
dataset = dss_file_obj.read_grid(pathname)
profile = {
'driver': 'GTiff',
'dtype': dataset.profile['opt_dtype'],
'width': dataset.width,
'height': dataset.height,
'count': 1,
'crs': crs_raster,
'transform': dataset.transform,
'nodata': -999
}
with rasterio.open(save_file, 'w', **profile) as dst:
dst.write(dataset.read(), 1)

But yeah the solution you presented would be much easier.

Also one of the reason I have noticed the crs.to_wkt fails is because the dss grid crs uses "metre" instead of "Meter" as units. if you still want to read the crs code from the dss grid one way that works is by replacing it from the string.

crs_raster = CRS.from_wkt(dataset.crs.replace('["metre",0.0]', '["Meter",1.0]'))

mxkpp commented

@KushalMBI glad there is more than one way to get around this. Thanks for sharing.