Skip to content

Raster reprojection in R

3 messages · Víctor Rodríguez Galiano, Michael Sumner

Original
Víctor Rodríguez Galiano
# 
Hello,

 I am trying to reproject a raster image from sinusoidal projection using
?projectRaster?. The size of the image is not very big (922 KB), but when
applying the reprojection I get this error message: ?Error: cannot allocate
vector of size 7.1 Gb?. Please see the code below:

HDFpath <- "C:/images/" # dir with images
setwd(HDFpath)            # set working directory
library(raster)
library(rgdal)
geocrs <- "+proj=longlat +ellps=WGS84 +datum=WGS84"
tile <- brick(?image.tif?)
tile_reproj <- projectRaster(tile, crs=geocrs)

Error: cannot allocate vector of size 7.1 Gb
class       : RasterStack 
dimensions  : 2400, 2400, 5760000, 1  (nrow, ncol, ncell, nlayers)
resolution  : 463.3127, 463.3127  (x, y)
extent      : -11119505, -10007555, 5559753, 6671703  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs 
names       : INCA.h08v03.Dormancy_median 
min values  :                          19 
max values  :                         540
Michael Sumner
# 
You should set a target raster with the extent and dimensions required.
There are inherent limits in reprojection and heuristics won't always work.
Generally using a target raster is much more efficient anyway.

But, this is an intensive remodeling of the data, delivered in a projection
for good reason (global equal area probably, and reasonably close to the L3
bins used for daily statistics.

You should find an alternative process IMO, are you trying to extract pixel
values or something else?

Cheers, Mike

On Sat., 14 Sep. 2019, 11:30 V?ctor Rodr?guez Galiano, <vrgaliano at gmail.com>
wrote:

  
  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

      
Hi again, I checked the actual extent and it makes more sense now as it is
a tiny region of the Aleutians. The problem in raster's heuristic is that
part of your projected raster overlaps the anti-meridan, and so the
extent-determiner unhelpfully expands to include the full extent of all
longitudes:

library(raster)
ex <- extent(-11119505, -10007555, 5559753, 6671703)
prj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+ +b=6371007.181 +units=m +no_defs"
projectExtent(raster(ex, crs = prj, res = 463), crs = "+init=epsg:4326")
class      : RasterLayer
dimensions : 2402, 2402, 5769604  (nrow, ncol, ncell)
resolution : 0.1498734, 0.004163854  (x, y)
extent     : -179.9968, 179.999, 49.99842, 60  (xmin, xmax, ymin, ymax)
crs        : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0

Illustrated here:
https://gist.github.com/mdsumner/012c4c8e36ffa42e9053621420924eab

So, you'd want something like

projectRaster(tile, raster(extent(-180, -135, 50, 60), crs = geocrs, res =
0.02)

but, that's still outrageously wasteful and unnecessary as mentioned before
- but depends on what you are doing.

(It's a curious way to store data, but I'm pretty sure the MODIS terra
community settled on a global raster in sinusoidal to match the L3 bins in
a way the marine community doesn't usually do).

Thanks for the example, that gives a good way to approach this question on
various forums where it comes up quite a lot and I never was really
motivated to pursue before.

Cheers, Mike.

        
On Sat, Sep 14, 2019 at 7:49 PM Michael Sumner <mdsumner at gmail.com> wrote: