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