Cannot use PROJ4 with gdal.Translate in Python

Question

I'm trying to convert MCD12Q1 HDF5 product to GeoTiff. I extracted metadata from the file, and from this user guide I found the projection information as

PROJ4: '+proj=sinu +a=6371007.181 +b=6371007.181 +units=m'

However, when I tried to use these variables in my Python script:

from osgeo import gdal
 
Dataset = gdal.Open("../data/MCD12Q1/MCD12Q1.A2001001.h19v04.006.2018142205330.hdf", gdal.GA_ReadOnly)
SubDatasets = Dataset.GetSubDatasets()
SubLayers = [SubLayer for SubLayer, _ in SubDatasets]
SubLayer = gdal.Open(SubLayers[0], gdal.GA_ReadOnly)
SLName = SubLayer.GetMetadata_Dict()['long_name']
print(SLName)
WestCoord = SubLayer.GetMetadata_Dict()["WESTBOUNDINGCOORDINATE"]
NorthCoord = SubLayer.GetMetadata_Dict()["NORTHBOUNDINGCOORDINATE"]
EastCoord = SubLayer.GetMetadata_Dict()["EASTBOUNDINGCOORDINATE"]
SouthCoord = SubLayer.GetMetadata_Dict()["SOUTHBOUNDINGCOORDINATE"]
PROJ4 = "+proj=sinu +a=6371007.181 +b=6371007.181 +units=m"
Command = f"-a_srs '{PROJ4}' -a_ullr {WestCoord} {NorthCoord} {EastCoord} {SouthCoord}"
print(Command)
Options = gdal.TranslateOptions(gdal.ParseCommandLine(Command))
gdal.Translate("../data/MCD12Q1/MCD12Q1.A2001001.h19v04.006.2018142205330_.tif", SubLayer, options=Options)
>>> Land_Cover_Type_1
>>> -a_srs '+proj=sinu +a=6371007.181 +b=6371007.181 +units=m' -a_ullr 13.054073 50.0 31.127441 40.0
>>> Traceback (most recent call last):
 File "D:\geo\script\hdf5_to_geotiff.py", line 21, in <module>
 gdal.Translate("../data/MCD12Q1/MCD12Q1.A2001001.h19v04.006.2018142205330_.tif", SubLayer, options=Options)
 File "C:\ProgramData\Miniconda3\envs\geo\lib\site-packages\osgeo\gdal.py", line 422, in Translate
 return TranslateInternal(destName, srcDS, opts, callback, callback_data)
 File "C:\ProgramData\Miniconda3\envs\geo\lib\site-packages\osgeo\gdal.py", line 3381, in TranslateInternal
 return _gdal.TranslateInternal(*args)
TypeError: in method 'TranslateInternal', argument 3 of type 'GDALTranslateOptions *'

In this script, EPSG code is used but I only have PROJ4 and WKT formats. How can I use PROJ4 or WKT definitions with gdal.Translate function?

• snowman2: Does this answer your question: gis.stackexchange.com/questions/323317/…
• snowman2: You can try: rioxarray.open_rasterio as well.
• mdsumner: use +R not a,b .

Answer 1

This is the final version of my script. I still cannot understand how it worked, though. I don't think I used the projected coordinates in the correct order.

from osgeo import gdal
from pyproj import Proj
Dataset = gdal.Open("../data/MCD12Q1/MCD12Q1.A2001001.h19v04.006.2018142205330.hdf", gdal.GA_ReadOnly)
SubDatasets = Dataset.GetSubDatasets()
SubLayers = [SubLayer for SubLayer, _ in SubDatasets]
SubLayer = gdal.Open(SubLayers[0], gdal.GA_ReadOnly)
SLName = SubLayer.GetMetadata_Dict()['long_name']
proj4 = "+proj=sinu +a=6371007.181 +b=6371007.181 +units=m"
crs_prj = Proj(proj4)
ulcLat = float(SubLayer.GetMetadata_Dict()["NORTHBOUNDINGCOORDINATE"])
ulcLon = float(SubLayer.GetMetadata_Dict()["WESTBOUNDINGCOORDINATE" ])
lrcLat = float(SubLayer.GetMetadata_Dict()["SOUTHBOUNDINGCOORDINATE"])
lrcLon = float(SubLayer.GetMetadata_Dict()["EASTBOUNDINGCOORDINATE" ])
lrcX, ulcY = crs_prj(lrcLon, ulcLat)
ulcX, lrcY = crs_prj(ulcLon, lrcLat)
print(ulcX, ulcY, lrcX, lrcY)
gdal.Translate("../data/MCD12Q1/MCD12Q1.A2001001.h19v04.006.2018142205330_.tif",
 SubLayer,
 **{"outputBounds": (ulcX, ulcY, lrcX, lrcY),
 "outputSRS" : proj4})

enter image description here

Answer 2

What you provide does not match the API (or it will not throw an error due TypeError e.g https://docs.python.org/3/library/exceptions.html#TypeError)

Instead of providing a third argument using gdal.TranslateOptions, you can directly provide the args according to the API doc https://gdal.org/python/osgeo.gdal-module.html#Translate

Keyword arguments are : options --- return of gdal.TranslateOptions(), string or array of strings other keywords arguments of gdal.TranslateOptions() If options is provided as a gdal.TranslateOptions() object, other keywords are ignored.

So, try

gdal.Translate("../data/MCD12Q1/MCD12Q1.A2001001.h19v04.006.2018142205330_.tif", SubLayer, options=Command)

PS: I've already documented the way to provide options for gdal.VectorTranslate and the logic is similar for gdal.Translate. It may help you. If curious, you can take a look at Issue to convert from PostgreSQL input to GPKG using Python GDAL API function gdal.VectorTranslate

• gokdumano: it took me a while but i think i solved my problem. Your explanation helped me to use gdal.Translate. However there was another problem, I realized that I needed to transform geodetic coordinates to projected coordinates. The script below somehow worked but