Raster 02: When Rasters Don‘t Line Up - Reproject Raster Data in R

#raster02 When Rasters Don't Line Up - Reproject Raster Data in R
library(pacman)
p_load(raster, rgdal)
wd <-"G:/Rdata/neon_data/NEONDSAirborneRemoteSensing/"
setwd(wd)

#import DTM
DTM_HARV <- raster(paste0(wd, "HARV/DTM/HARV_dtmCrop.tif"))

#import DTM hillshade
DTM_hill_HARV <- raster(paste0(wd,"HARV/DTM/HARV_DTMhill_WGS84.tif"))

plot(DTM_hill_HARV, 
     col=grey(1:100/100),
     legend=FALSE,
     main="DTM Hillshade\n NEON Harvard Forest Field Site")

#overlay the DTM on the top of the hillshade
plot(DTM_HARV,
     col=terrain.colors(100),
     alpha = 1,
     legend=FALSE,
     main="Digital Terrain Model\nNEON Harvard Forest Field Site")

#发现无法叠加上去，提示内存不足，我把add = true给删去才可以绘图
#查看是否是坐标参考系不一致导致的无法叠加
crs(DTM_HARV)
# CRS arguments:
#         +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
crs(DTM_hill_HARV)
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs 

#果真是坐标参考系不一致：一个是投影坐标系，一个是地理坐标系
#请记住，仅当您首先为要重新投影的栅格对象定义了CRS时，重新投影才起作用。
#reproject to UTM

DTM_hill_HARV_UTM <- projectRaster(DTM_hill_HARV,crs=crs(DTM_HARV))
crs(DTM_HARV)
# CRS arguments:
#         +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
crs(DTM_hill_HARV)
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs 
crs(DTM_hill_HARV_UTM)
# CRS arguments:
#         +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 

extent(DTM_hill_HARV)
# class      : Extent 
# xmin       : -72.18192 
# xmax       : -72.16061 
# ymin       : 42.52941 
# ymax       : 42.54234

extent(DTM_hill_HARV_UTM)
# class      : Extent 
# xmin       : 731397.3 
# xmax       : 733205.3 
# ymin       : 4712403 
# ymax       : 4713907

extent(DTM_HARV)
# class      : Extent 
# xmin       : 731453 
# xmax       : 733150 
# ymin       : 4712471 
# ymax       : 4713838 

#目前经过重投影，空间参考系是一致了，但是范围还是不一致
#下面看下分辨率
res(DTM_HARV)
# [1] 1 1
res(DTM_hill_HARV)  #很明显是单位是经纬度
# [1] 1.22e-05 8.99e-06
res(DTM_hill_HARV_UTM)  #很奇怪，像素居然不是正方形的，有意思了
# [1] 1.000 0.998
#force our newly reprojected raster to be 1m x 1m resolution by add res=1
DTM_hill_HARV_UTM <- projectRaster(DTM_hill_HARV,crs=crs(DTM_HARV), res=1)
res(DTM_hill_HARV_UTM)
# [1] 1 1

#plot newly reprojected hillshade
plot(DTM_hill_HARV_UTM,
     col=grey(1:100/100),
     legend = F,
     main = "DTM with Hillshade\n NEON Harvard Forset Field Site")
plot(DTM_HARV,
     col=terrain.colors(100),
     alpha=0.4,
     add=T,
     legend=F)