0

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?

3

2 Answers 2

0

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

1
  • 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
    Commented Jul 1, 2021 at 3:02
0

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

Your Answer

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

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