1

I obtained monthly TRMM data (TMPA/3B43) from here. It is in .HDF format for each month. I need it for other algorithm in R that uses netCDF file on yearly basis, so I wrote a code that reads .HDF file, converts it to .tif file (with the help of James' answer), then reads it as raster and converts it to matrix.

The produced raster file (rast) doesn't have a geographical reference. What I understood so far is I have to assign the crs from the southwest corner to (-180,-50) with an interval 0.25 deg in both direction.

Latitude should be from the left: -50 to 49.75 by 0.25 will result in 400 columns Longitude should be from the bottom -180 to 179.75 by 0.25 will result in 1440 rows

My Question: how to assign this type of CRS to the raster in R? Below are the information of .HDF, .tiff, and matrix.

library(gdalUtils)
library(raster)

#use R wrapper for the 'gdalinfo' function to obtain metadata
> gdalinfo("3B43.19980101.7.HDF")
[1] "Driver: HDF4/Hierarchical Data Format Release 4"                                                  
[2] "Files: 3B43.19980101.7.HDF"                                    
[3] "Size is 512, 512"                                                                                 
[4] "Metadata:"                                                                                        
[5] "  FileHeader=AlgorithmID=3B43;"                                                                   
[6] "AlgorithmVersion=3B43_7.0;"                                                                       
[7] "FileName=3B43.19980101.7.HDF;"                                                                    
[8] "GenerationDateTime=2012-05-14T16:29:45.000Z;"                                                     
[9] "StartGranuleDateTime=1998-01-01T00:00:00.000Z;"                                                   
[10] "StopGranuleDateTime=1998-01-31T23:59:59.999Z;"                                                    
[11] "GranuleNumber=;"                                                                                  
[12] "NumberOfSwaths=0;"                                                                                
[13] "NumberOfGrids=1;"                                                                                 
[14] "GranuleStart=;"                                                                                   
[15] "TimeInterval=MONTH;"                                                                              
[16] "ProcessingSystem=PPS;"                                                                            
[17] "ProductVersion=7;"                                                                                
[18] "MissingData=;"                                                                                    
[19] ""                                                                                                 
[20] "  FileInfo=DataFormatVersion=m;"                                                                  
[21] "TKCodeBuildVersion=1;"                                                                            
[22] "MetadataVersion=m;"                                                                               
[23] "FormatPackage=HDF Version 4.2 Release 4, January 25, 2009;"                                       
[24] "BlueprintFilename=TRMM.V7.3B43.blueprint.xml;"                                                    
[25] "BlueprintVersion=BV_13;"                                                                          
[26] "TKIOVersion=1.6;"                                                                                 
[27] "MetadataStyle=PVL;"                                                                               
[28] "EndianType=LITTLE_ENDIAN;"                                                                        
[29] ""                                                                                                 
[30] "  GridHeader=BinMethod=ARITHMETIC_MEAN;"                                                          
[31] "Registration=CENTER;"                                                                             
[32] "LatitudeResolution=0.25;"                                                                         
[33] "LongitudeResolution=0.25;"                                                                        
[34] "NorthBoundingCoordinate=50;"                                                                      
[35] "SouthBoundingCoordinate=-50;"                                                                     
[36] "EastBoundingCoordinate=180;"                                                                      
[37] "WestBoundingCoordinate=-180;"                                                                     
[38] "Origin=SOUTHWEST;"                                                                                
[39] ""                                                                                                 
[40] "Subdatasets:"                                                                                     
[41] "  SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:\"/BE70E3A0F5C15471/TRMM/3B43.19980101.7.HDF\":0"
[42] "  SUBDATASET_1_DESC=[1440x400] precipitation (32-bit floating-point)"                             
[43] "  SUBDATASET_2_NAME=HDF4_SDS:UNKNOWN:\"/BE70E3A0F5C15471/TRMM/3B43.19980101.7.HDF\":1"
[44] "  SUBDATASET_2_DESC=[1440x400] relativeError (32-bit floating-point)"                             
[45] "  SUBDATASET_3_NAME=HDF4_SDS:UNKNOWN:\"/BE70E3A0F5C15471/TRMM/3B43.19980101.7.HDF\":2"
[46] "  SUBDATASET_3_DESC=[1440x400] gaugeRelativeWeighting (8-bit integer)"                            
[47] "Corner Coordinates:"                                                                              
[48] "Upper Left  (    0.0,    0.0)"                                                                    
[49] "Lower Left  (    0.0,  512.0)"                                                                    
[50] "Upper Right (  512.0,    0.0)"                                                                    
[51] "Lower Right (  512.0,  512.0)"                                                                    
[52] "Center      (  256.0,  256.0)"  

#Returns HDF subdataset names for standardized files.
> sds <- get_subdatasets("3B43.19980101.7.HDF")

# precipiction is the first subset, convert to tif
> gdal_translate(sds[1], dst_dataset = "NPP2000.tif")


       
# Load the new .tif
       
> rast <- raster("NPP2000.tif") 

> rast
class      : RasterLayer 
dimensions : 1440, 400, 576000  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 0, 400, 0, 1440  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : /media/ahmed/Volume/ITCZ_algo_dev/ITCZ_Algorithm/NPP2000.tif 
names      : NPP2000 

> med<-raster::as.matrix(rast)
> dim(med)
[1] 1440  400
> plot(rast)

The plotted raster layer

2

You can use terra like this

library(terra)
# terra version 1.2.5
x <- rast("3B43.19980101.7.HDF")
#Warning message:
#[rast] unknown extent
 
x
#class       : SpatRaster 
#dimensions  : 1440, 400, 3  (nrow, ncol, nlyr)
#resolution  : 0.0025, 0.0006944444  (x, y)
#extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#sources     : 3B43.19980101.7.HDF:0  
#              3B43.19980101.7.HDF:1  
#              3B43.19980101.7.HDF:2  
#names       : 0, 1, 2 

Now fix the data

y <- t(x)
z <- flip(y, "v")
ext(z) <- c(-180, 180, -50, 50)
z

#class       : SpatRaster 
#dimensions  : 400, 1440, 3  (nrow, ncol, nlyr)
#resolution  : 0.25, 0.25  (x, y)
#extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       :           0,           1,           2 
#min values  :  0.00000000,  0.00144434,  0.00000000 
#max values  :   1.5933868,   0.5838504, 100.0000000 

plot(z)

enter image description here

If you only care about one of the layers you can do, for example to get the second one,

x <- rast("3B43.19980101.7.HDF", 2)
z <- flip(t(x), "v")
ext(z) <- c(-180, 180, -50, 50)
3
  • You actually saved my life! Many thanks! Very much appreciated!
    – ahmathelte
    Apr 28 at 19:37
  • I have another Q: I try to install your nice libraray by: install.packages("terra", configure.args = "--with-gdal-config=/home/ahmed/miniconda2/pkgs/libgdal-3.1.2-hb2a6f5f_2/bin/gdal-config"). It stops at: checking for gdal.h... no configure: error: gdal.h not found in given locations. ERROR: configuration failed for package ‘terra’ I have miniconda installed.
    – ahmathelte
    Apr 28 at 20:27
  • 1
    I assume you need to install GDAL first. Here are some instructions for installation: github.com/rspatial/terra Apr 28 at 20:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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