Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Why the values of my raster map change when I project it to a new CRS (projectRaster function in R)?

I need to project a map in latlong to an azimuthal equidistant projection.

map_proj<-projectRaster(map, crs="+proj=aeqd +lon_0=48+lat_0=18")

In my original map I have these values

class       : RasterLayer 
dimensions  : 1713, 2185, 3742905  (nrow, ncol, ncell)
resolution  : 0.008335028, 0.00833354  (x, y)
extent      : 39.06984, 57.28187, -25.93814, -11.66279  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /Users/Oritteropus/Desktop/Progetti/maps/mada_rast2 
names       : mada_rast2 
values      : 0, 255  (min, max)

unique(map)
 [1]  0  1  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 24 25 26 27 28
[24] 29 30 31 32 33 34 35 36 37 38 39 40

In my projected map I now have the following values

class       : RasterLayer 
dimensions  : 1990, 2254, 4485460  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : 4026994, 6280994, -3379165, -1389165  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aeqd +lon_0=0+lat_0=0 +ellps=WGS84 
data source : in memory
names       : mada_rast2 
values      : -0.1498806, 40  (min, max)

head(unique(map_proj))

[1] -1.498806e-01 -8.050509e-02  0.000000e+00  1.164434e-05  3.575309e-05
[6]  5.233506e-05

I need to keep the same values.. Can somebody explain me what happens or where I'm wrong? Thanks

like image 793
Oritteropus Avatar asked Mar 26 '13 10:03

Oritteropus


People also ask

How do you change the CRS of a raster?

We can use the projectRaster() function to reproject a raster into a new CRS. Keep in mind that reprojection only works when you first have a defined CRS for the raster object that you want to reproject. It cannot be used if no CRS is defined.

How do I project rasters in R?

There are two approaches you can follow to project the values of a Raster object. 1) Provide a crs argument, and, optionally, a res argument, but do not provide a to argument. 2) Create a template Raster with the CRS you want to project to.


1 Answers

It is because re-projecting a raster requires a certain amount of warping of the cells and by necessity interpolation between the values. If you want to keep the same values you will have to change the method of interpolation. Use nearest neighbour instead...

map_proj<-projectRaster(map, crs="+proj=aeqd +lon_0=48+lat_0=18" , method = "ngb" )

An example with data from the help page for raster:::projectRaster:

r <- raster(xmn=-110, xmx=-90, ymn=40, ymx=60, ncols=40, nrows=40)
r <- setValues(r, 1:ncell(r))
newproj <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"

# Default reproject uses bilinear interpolation
r_bl <- projectRaster( r , crs = newproj )
r_nn <- projectRaster( r , crs = newproj , method = "ngb" )

range(values(r) , na.rm = T) 
#[1]    1 1600
range(values(r_bl) , na.rm = T )
#[1]   -7.671638 1612.968493
range(values(r_nn) , na.rm = T )
#[1]    1 1600

Of course you are still projecting values to new cells so you will probably get some NA's. Briefly, bilinear interpolation makes sense for continuous variables, whilst nearest neighbour makes sense for categorical variables.

like image 181
Simon O'Hanlon Avatar answered Oct 25 '22 23:10

Simon O'Hanlon