[R-sig-Geo] Raster reprojection in R

[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 using gmail.com> 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, <
> vrgaliano using gmail.com> 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
>> R-sig-Geo using r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>

-- 
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsumner using gmail.com

	[[alternative HTML version deleted]]