0

I am interested by MODIS product. I've converted HDF file to TIFF with R and gdalUtils. I used a method similar to this : https://stackoverflow.com/questions/36772341/reading-hdf-files-into-r-and-converting-them-to-geotiff-rasters

But I found that the product is in sinusoidal projection, and I want to convert it to WGS84. I have tested with R, but files is maybe too big (too many cells because it is a high resolution, more details of the product here https://lpdaac.usgs.gov/products/mod17a3hgfv006/)

I use the following code, where I'm loading first the convert tiff file (and found TIFF fill for the specific year, here 2000). Then I created an empty raster with world extent and define a projection. I finally use mosaic_rasters, to combine all the different TIFF files together.

listtiff <- list.files(, pattern = "\\.tif$", full.names = TRUE) 
  
subfileyear <- listtiff[grepl(paste0(2000), listtiff)] #find file with pattern "2000"
  
e <- raster::extent(-180.0000, 180,  -90.0000  , 90.0000)
template <- raster::raster(e)
#proj4string(template) <- CRS("+init=epsg:4326")
proj4string(template) <- CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
writeRaster(template, file="NPP2000.tif", format="GTiff", overwrite = TRUE)
mosaic_rasters(gdalfile=subfileyear,dst_dataset="NPP2000.tif",of="GTiff" , gdalwarp_index = 1)

If I'm using gdal info I obtain this, but now I do not know how to convert it to a wgs84??

I tried to use projectRaster but it did not work well, I think it's due to the large number of pixels. Do you know other option for conversion from ellipsoide to WGS84?

> gdalinfo("NPP2000.tif")
 [1] "Driver: GTiff/GeoTIFF"                                                     
 [2] "Files: NPP2000.tif"                                                        
 [3] "Size is 86400, 36000"                                                      
 [4] "Coordinate System is:"                                                     
 [5] "PROJCRS[\"unnamed\","                                                      
 [6] "    BASEGEOGCRS[\"Unknown datum based upon the custom spheroid\","         
 [7] "        DATUM[\"Not_specified_based_on_custom_spheroid\"," 
 [8] "            ELLIPSOID[\"Custom spheroid\",6371007.181,0,"                  
 [9] "                LENGTHUNIT[\"metre\",1,"                                   
[10] "                    ID[\"EPSG\",9001]]]],"                                 
[11] "        PRIMEM[\"Greenwich\",0,"                                           
[12] "            ANGLEUNIT[\"degree\",0.0174532925199433,"                      
[13] "                ID[\"EPSG\",9122]]]],"                                     
[14] "    CONVERSION[\"Sinusoidal\","                                            
[15] "        METHOD[\"Sinusoidal\"],"                                           
[16] "        PARAMETER[\"Longitude of natural origin\",0,"                      
[17] "            ANGLEUNIT[\"degree\",0.0174532925199433],"                     
[18] "            ID[\"EPSG\",8802]],"                                           
[19] "        PARAMETER[\"False easting\",0,"                                    
[20] "            LENGTHUNIT[\"metre\",1],"                                      
[21] "            ID[\"EPSG\",8806]],"                                           
[22] "        PARAMETER[\"False northing\",0,"                                   
[23] "            LENGTHUNIT[\"metre\",1],"                                      
[24] "            ID[\"EPSG\",8807]]],"                                          
[25] "    CS[Cartesian,2],"                                                      
[26] "        AXIS[\"easting\",east,"                                            
[27] "            ORDER[1],"                                                     
[28] "            LENGTHUNIT[\"metre\",1,"                                       
[29] "                ID[\"EPSG\",9001]]],"                                      
[30] "        AXIS[\"northing\",north,"                                          
[31] "            ORDER[2],"                                                     
[32] "            LENGTHUNIT[\"metre\",1,"                                       
[33] "                ID[\"EPSG\",9001]]]]"                                      
[34] "Data axis to CRS axis mapping: 1,2"                                        
[35] "Origin = (-20015109.353999998420477,10007554.676999999210238)"             
[36] "Pixel Size = (463.312716527770533,-463.312716527777638)"                   
[37] "Metadata:"                                                                 
[38] "  AREA_OR_POINT=Area"                                                      
[39] "Image Structure Metadata:"                                                 
[40] "  INTERLEAVE=BAND"                                                         
[41] "Corner Coordinates:"                                                       
[42] "Upper Left  (-20015109.354,10007554.677) (124d17'29.21\"W, 90d 0' 0.00\"N)"
[43] "Lower Left  (-20015109.354,-6671703.118) (  0d 0' 0.00\"E, 60d 0' 0.00\"S)"
[44] "Upper Right (20015109.354,10007554.677) (124d15' 5.75\"E, 90d 0' 0.00\"N)" 
[45] "Lower Right (20015109.354,-6671703.118) (  0d 0' 0.00\"W, 60d 0' 0.00\"S)" 
[46] "Center      (      -0.000, 1667925.780) (  0d 0' 0.00\"W, 15d 0' 0.00\"N)" 
[47] "Band 1 Block=86400x1 Type=Int16, ColorInterp=Gray"                         
[48] "  NoData Value=32767"                                                      
[49] "  Offset: 0,   Scale:0.0001"
2
  • 2
    What does "it did not work well" mean? Did you get some errors? Please edit your question and show us what went wrong. But yes, that is a very large raster.
    – Spacedman
    Feb 4 at 19:46
  • It takes a lot of times, and R crash (I obtain a message R session aborted). I tried with only one year, and when it works (after restarting R session after loading the raster for example) I have a similar issue to this person, where I have only data in the "northern" hemisphere : stackoverflow.com/questions/47800587/… The problem is that I do not want to aggregate the data to a lower resolution
    – Vale
    Feb 4 at 22:31
1

Just taking a shot here.. I've worked with big raster datasets a lot recently, and while I mostly use R, for converting to a new projection I found that QGIS is much more efficient. So my advice would be to try to open your raster in QGIS and use "Raster"->"Projections"->"Warp (Reproject)". You can then save your re-projected raster, load it into R and keep going from there..

1

Using the answer of @Lena, I maybe found a R solution with the gdalUtils

You can directly reproject it with the following code, where you specify your raster file, the destination file (where outdire is the path where I want to put the results), with the projection of the MODIS product and the desire projection (here WGS84). Faster if you have mutlpile rasters to convert, and faster than raster package

gdalwarp("NPP2000.tif",dstfile=file.path(outdir,"NPP2000.tif"),s_srs = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs",
         t_srs='+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0',output_Raster=TRUE, overwrite=TRUE,verbose=TRUE) 
0

For Python users, I have written a short article and Gist explaining how to perform this translation in GDAL and Rasterio.

from osgeo import gdal
  
in_file = f"./modis.hdf" # raw MODIS HDF in sinusoid projection
out_file = f"./modis.tif"
    
# open dataset
dataset = gdal.Open(in_file,gdal.GA_ReadOnly)
subdataset =  gdal.Open(dataset.GetSubDatasets()[0][0], gdal.GA_ReadOnly)

# gdalwarp
kwargs = {'format': 'GTiff', 'dstSRS': 'EPSG:4326'}
ds = gdal.Warp(destNameOrDestDS=out_file,srcDSOrSrcDSTab=subdataset, **kwargs)
del ds

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.