2

I am trying to open, perform a calculation then save to a new tif file using RasterIO without changing the coordinate system.

My input tif is in GDA 1994 MGA Zone 54 and has a wkt string of

 PROJCS["GDA_1994_MGA_Zone_54",GEOGCS["GCS_GDA_1994",DATUM["GDA_1994",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",141],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

RasterIO, for the same tif, reads and writes the coordinate system as UTM Zone 54, Southern Hemisphere and which has a wkt string of

PROJCS["UTM Zone 54, Southern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",141],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]

Both these coordinate system have the same Proj4 string. It seems that RasterIO creates its crs.wkt string from the proj4 string, not from the one provided by the input tif.

Is there a way using RasterIO to read the wkt string from the tif and use it when saving to a new file?

Here's my code so far:

smin=0
smax=255
with rasterio.open('input_94mga54.tif') as src:
    band = src.read(1,masked=True)
    rescaled = (band - band.min()) * ( smax - smin) / (band.max() - band.min()) + smin
    with rasterio.open('rescaled_94mga54.tif', 'w', **src.meta) as dst:
        dst.write_band(1, rescaled)

1 Answer 1

1

Since posting this question, I have found a partial answer.

You can assign a WKT string of the coordinate system using

from osgeo import gdal
out_meta['crs'] = gdal.Open('input_94mga54.tif').GetProjectionRef()

This saves the correct WKT to the file (confirmed by checking using GDAL methods) however as soon as you read the new file and look at the WKT using src.crs.wkt it displays the incorrect WKT of UTM Zone 54, Southern Hemisphere.

Not sure what impact has on reprojections.


Just another update - This has now been fixed in rasterio v1.1.14 and above.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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