Raster reprojection in R

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Raster reprojection in R

Victor Rodriguez-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

> tile
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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Raster reprojection in R

Michael Sumner-2
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, <[hidden email]>
wrote:

> 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
>
> > tile
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Raster reprojection in R

Michael Sumner-2
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 <[hidden email]> wrote:

> 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, <
> [hidden email]> wrote:
>
>> 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
>>
>> > tile
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>

--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: [hidden email]

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo