Geographic Information Systems Asked by Jorge V on August 24, 2021
I downloaded a complete month of daily images of MODIS MOD07_L2
product (Atmospheric profile). I’ve worked with *.hdf
files for really long (using both R
and tools such as HEG
), nevertheless, I am now stuck with MOD07 product. I need to work with Dew-point Temperature at different pressure levels, but I don’t understand how the data is actually saved in the *.hdf
layer (which is Retrieved_Moisture_Profile
). I tried using HEG
software with all possible options (for output file types), but none of them worked. For every iteration I tried in HEG I got the "Could not find TimeofDay attribute, CalendarDate, Sensorname, pointingangle (and many more) attribute in input hdf file". So I believe HEG
is not capable to work with this type of data. The result for every image was a black square with just one pixel value.
I opened the *.hdf
file directly in QGIS
(version 3.10.4-A Coruña) and, more specifically, I opened the Retrieved_Moisture_Profile
layer. When I studied the histogram of the latter image it is quite clear that there are 20 different Dew-point temperature profiles (each of one of them is for an specific pressure level), as you can see in the following figure.
I was thinking that maybe I can just save the *.hdf
to *.nc
, but I tried with rhdf5
and the latter package is depricated for R version 3.6.3
. I also tried opening the *.hdf
file with gdalUtils
, which actually works, but I couldn’t find how to ‘divide’ each layer of the raster or even save them separately.
So, what I need is the following (in R
).
*.hdf
fileRetrieved_Moisture_Profile
layer within the *.hdf
fileGeoTIFF
file.FYI: I’m running R version 3.6.3
in Windows 10 64 bit machine.
Thanks in advance for your help!
Sincerely, Jorge.
To anyone who is struggling with MOD07_L2
product and was wondering how I solved it, I followed this answer from @fdetsch to Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal, but modified a little bit. So basically what I did is the following.
library(raster)
library(rgdal)
library(gdalUtils)
# hdf files directory
main_dir <- file.path("D:/path_to_your_hdf_files...")
# List all *.hdf files within main_dir
hdf_files <- list.files(path = main_dir, pattern = "*.hdf", full.names = T)
# Get the info of the first hdf within the hdf_files list
hdf_info <- GDALinfo(hdf_files[1], returnScaleOffset = FALSE)
# Obtain the metadata of the latter hdf file
meta <- attr(hdf_info, "mdata")
# Name the extent name attributes to search
crd_str <- paste0(c("WEST", "EAST", "SOUTH", "NORTH"), "BOUNDINGCOORDINATE")
# Function to search for the attributes
crd_id <- sapply(crd_str, function(i) grep(i, meta))
# Extract the coordinates of each bounding coordinate
crd <- meta[crd_id]
crd <- as.numeric(sapply(strsplit(crd, "="), "[[", 2))
# Define the coordinates as a raster::extent to later add to the raster
ext <- extent(crd)
# Open the first hdf file
one <- getSds(hdf_files[1])
# Open band 16 (which is Retrieved Moisture Profile) as a raster brick
b16 <- brick(readGDAL(one$SDS4gdal[16], as.is = TRUE))
# Set crd as the extent of b16
extent(b16) <- crd
# Set EPSG:4326 as the projection of b16
projection(b16) <- "+init=epsg:4326"
The result will be a RasterBrick
object with nlayer = 20
and each layer (which is an specific pressure level) will have the dew point temperature values (however you must add the offset and scale factor).
Answered by Jorge V on August 24, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP