12

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

J. Win.
  • 6,662
  • 7
  • 34
  • 52
Oritteropus
  • 531
  • 2
  • 6
  • 18

1 Answers1

16

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.

Simon O'Hanlon
  • 58,647
  • 14
  • 142
  • 184