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?

Does this answer your question: https://gis.stackexchange.com/questions/323317/converting-netcdf-dataset-array-to-geotiff-using-rasterio-python/329141#329141 — snowman2, Jun 30 '21 at 21:26

ThomasG77 · Answer 1 · 2021-06-30T22:09:06.380

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

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 — gokdumano, Jul 01 '21 at 03:02

score 0 · Answer 2 · edited Dec 15 '22 at 13:18

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})

Cannot use PROJ4 with gdal.Translate in Python

2 Answers2