Converting sinusoidal to WGS - MODIS TIFF file too big with R

Question

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"

Answer 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..

Answer 2

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) 

Answer 3

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