TransWikia.com

GRASS GIS: r.in.gdal does not find raster data in .hdf file (modis data)

Geographic Information Systems Asked by F.H. on December 9, 2020

I’m having trouble in importing MODIS .hdf data in GRASS GIS. When I use r.in.gdal it returns me this error:

ERROR: No raster bands found in </Volumes/Maxtor/TESI_satellite_imagery/IMMAGINI_DI_PARTENZA/serie_storica_estate_dic2015-feb2020/modis/MOD02QKM.A2019044.1415.061.2019045012150.hdf>.

I’ve tried to uninstall and reinstall gdal on my computer (S.O. MacOS Catalina), and rebuild all the bindings with python, but I didn’t manage to solve the issue. I used Homebrew and followed its instructions to do it. I’ve also tried to use conda previously, but it has issues while "solving the environment" and cannot complete the procedure. What am I doing wrong?

I’ve also tried the i.modis.import, but I get this error

ERROR: pymodis library is not installed.

I can’t understand where this comes from, since if I run pip install pymodis or pip3 install pymodis on the terminal it says:

Requirement already satisfied: pymodis in ./opt/anaconda3/lib/python3.8/site-packages (2.0.9)
Requirement already satisfied: future in ./opt/anaconda3/lib/python3.8/site-packages (from pymodis) (0.18.2)
Requirement already satisfied: requests in ./opt/anaconda3/lib/python3.8/site-packages (from pymodis) (2.24.0)
Requirement already satisfied: GDAL in ./opt/anaconda3/lib/python3.8/site-packages (from pymodis) (3.1.3)
Requirement already satisfied: numpy in ./opt/anaconda3/lib/python3.8/site-packages (from pymodis) (1.18.5)
Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in ./opt/anaconda3/lib/python3.8/site-packages (from requests->pymodis) (1.25.9)
Requirement already satisfied: certifi>=2017.4.17 in ./opt/anaconda3/lib/python3.8/site-packages (from requests->pymodis) (2020.6.20)
Requirement already satisfied: idna<3,>=2.5 in ./opt/anaconda3/lib/python3.8/site-packages (from requests->pymodis) (2.10)
Requirement already satisfied: chardet<4,>=3.0.2 in ./opt/anaconda3/lib/python3.8/site-packages (from requests->pymodis) (3.0.4)

How can I solve this problem?

2 Answers

What works for the installation of gdal:

  1. install Anaconda
  2. created a separate environment: conda create -n gdal python=3.8
  3. activate the environment: conda activate gdal
  4. install inside the environment gdal: conda install -c conda-forge gdal
  5. exit the environment once finished using it: conda deactivate

The easiest way to use hdf files: install the HEG conversion tool by NASA. It allows both to convert the file into several formats and several reprojection options. Easiest way to access the metadata: install the HDFview tool by NASA. This worked better than gdal_translate, because the latter couldn't handle correctly some reprojection (it generated some weird artefacts, like horizontal stripes all over the tile), while the HEG tool managed to reproject all of my scenes correctly. I suggest to use the geographic reprojection, because using the UTM directly didn't gave a realistic result. I used the geographic and did a second rerojection within GRASS GIS.

I hope it helps anybody who'll face the same problems I had.

Correct answer by F.H. on December 9, 2020

Just to summarize from the comments:

gdalinfo shows the different HDF datasets as

Subdatasets:
  SUBDATASET_1_NAME=HDF4_EOS:EOS_SWATH:"MOD02QKM.A2019044.1415.061.2019045012150.hdf":MODIS_SWATH_Type_L1B:EV_250_RefSB
  SUBDATASET_1_DESC=[2x8120x5416] EV_250_RefSB MODIS_SWATH_Type_L1B (16-bit unsigned integer)
  SUBDATASET_2_NAME=HDF4_EOS:EOS_SWATH:"MOD02QKM.A2019044.1415.061.2019045012150.hdf":MODIS_SWATH_Type_L1B:EV_250_RefSB_Uncert_Indexes
  SUBDATASET_2_DESC=[2x8120x5416] EV_250_RefSB_Uncert_Indexes MODIS_SWATH_Type_L1B (8-bit unsigned integer)
  SUBDATASET_3_NAME=HDF4_SDS:UNKNOWN:"MOD02QKM.A2019044.1415.061.2019045012150.hdf":0
  SUBDATASET_3_DESC=[2030x1354] Latitude (32-bit floating-point)
  SUBDATASET_4_NAME=HDF4_SDS:UNKNOWN:"MOD02QKM.A2019044.1415.061.2019045012150.hdf":1
.....

Then you can import a single dataset into an unprojected XY Location with the command:

GRASS 7.8.4 (XY):~/work/tmp > r.in.gdal input=HDF4_SDS:UNKNOWN:"MOD02QKM.A2019044.1415.061.2019045012150.hdf":1 output=sds1
No projection information available
Importing raster map <sds1>...
 100%

If you want to reproject the MODIS datasets in advance then these steps might help: First use gdal_translate to create a VRT (should automatically contain the GCP's). Then gdalwarp to the projection that you want, using the "thin plate spline" option. And finally import into a GRASS Location defined by that same projection. i.e.

gdal_translate -of VRT HDF4_EOS:EOS_SWATH:"MOD02QKM.A2019044.1415.061.2019045012150.hdf":MODIS_SWATH_Type_L1B:EV_250_RefSB eos.vrt
gdalwarp -t_srs EPSG:4326 -tps eos.vrt eos.tif
# Switch to a WGS84 Location, and import
g.mapset location=WGS84 map=PERMANENT
r.in.gdal input=eos.tif output=eos
Importing 2 raster bands...
Importing raster map <eos.1>...
 100%
Importing raster map <eos.2>...
 100%
r.info eos.1
 +----------------------------------------------------------------------------+
 | Map:      eos.1                          Date: Tue Nov  3 20:52:04 2020    |
 | Mapset:   PERMANENT                      Login of Creator: micha           |
 | Location: WGS84                                                            |
 | DataBase: /home/micha/GIS/grass                                            |
 | Title:                                                                     |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    CELL                                                       |
 |   Rows:         7278                                                       |
 |   Columns:      9389                                                       |
 |   Total Cells:  68333142                                                   |
 |        Projection: Latitude-Longitude                                      |
 |            N: 12:08:47.681048S    S: 33:12:58.502358S   Res: 0:00:10.42193 |
 |            E: 47:24:22.942618W    W: 74:35:14.462442W   Res: 0:00:10.42193 |
 |   Range of data:    min = 367  max = 65534                                 |
 |                                                                            |
 |   Data Description:                                                        |
 |    generated by r.in.gdal                                                  |
 |                                                                            |
 |   Comments:                                                                |
 |    r.in.gdal input="eos.tif" output="eos" memory=300 offset=0 num_digit   |
 |    s=0                                                                     |
 |                                                                            |
 +----------------------------------------------------------------------------+

g.mapset location=WGS84 map=PERMANENT

Answered by Micha on December 9, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP