In QGIS with pyQgis I would like to determine if input raster has units in degrees or metres. Here I learned I can get projection from raster using GDAL:
inRaster = gdal.Open(rasterInPath,GA_ReadOnly)
projRef = inRaster.GetProjection()
print(projRef)
what's gives me string for raster with EPSG 4326:
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"]]
and for raster with EPSG 3857:
PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]
I guess I can do:
isInMetres = "PROJCS" in projRef
But is there any other solution? I don't consider this very readable and not sure if I can cover all corner cases this way (e.g. missing projRef). Dont need to be necessary with GDAL.
Using PyQgis API you can use the mapUnits method of QgsCoordinateReferenceSystem objects. It returns a QgsUnitTypes.DistanceUnit value, so you can use it like that for example :
# lyr is a QgsRasterLayer
crs = lyr.crs()
# crs is a QgsCoordinateReferenceSystem
unit = crs.mapUnits()
print('Unit is {}'.format(QgsUnitTypes.toString(unit)))
# 'Unit is mètres' (i'm using french language)
Reusing your code, you can write:
isInMetres = crs.mapUnits() == 0