R语言空间数据分析学习笔记——8栅格数据预处理

下载Landsat L2A数据:USGS EarthExplorer
现在USGS都可以之下下载大气校正后的产品了，以前还得[申请]（之前发的申请教程）(https://blog.csdn.net/qq_41718357/article/details/83903307)。
加载合肥市最新行政边界图

#栅格数据预处理
library(pacman)
p_load(sp,raster,rgeos,rgdal,rasterVis,raster,fs,sf,tidyverse,tmap)

#read data
hfbouds <-st_read(
  "G:\\Rdata\\Anhui\\hefei.shp"
)%>%
  st_transform(.,'+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs')
qtm(hfbouds,fill = "NAME")+
  tm_text("NAME")

下载DEM数据
裁剪：

hflandsat <- listlandsat %>%
  crop(.,hfbouds)%>%
  mask(.,hfbouds)
plot(hflandsat)

1
2
3
4
desert
bw
unicorn
使用imhof1风格添加水体

#栅格数据预处理
library(pacman)
p_load(sp,raster,rgeos,rgdal,rasterVis,raster,fs,sf,tidyverse,tmap,ggplot2,dplR,rayshader,plotly)

#read shp data
hfbouds <-st_read(
  "G:\\Rdata\\Anhui\\hefei.shp"
)%>%
  st_transform(.,'+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs')
crs(hfbouds)
qtm(hfbouds,fill = "NAME")+
  tm_text("NAME")

hfbouds1 <-st_transform(hfbouds,'+proj=longlat +datum=WGS84 +no_defs')


#read Landsat8 data
listlandsat<-dir("G:\\Rdata\\Anhui\\LC08_L2SP_121038_20201024_20201105_02_T1",
                 pattern="[B234567].TIF",full.names=T) %>% stack()#TIF一定要大写
Hf543 <-stack(listlandsat$LC08_L2SP_121038_20201024_20201105_02_T1_SR_B4,
              listlandsat$LC08_L2SP_121038_20201024_20201105_02_T1_SR_B3,
              listlandsat$LC08_L2SP_121038_20201024_20201105_02_T1_SR_B2)%>%
  crop(.,hfbouds)%>%
  mask(.,hfbouds)
plotRGB(Hf543,stretch='lin')


#裁剪
hflandsat <- listlandsat %>%
  crop(.,hfbouds)%>%
  mask(.,hfbouds)
plot(hflandsat)
# writeRaster(hflandsat,filename ="G:/Rdata/Anhui/LC08_L2SP_121038_20201024_20201105_02_T1/hfLandsat8.tif")


#read DEM data
hfDEM<- raster("G:\\Rdata\\Anhui\\DEM\\hfDEM.tif")%>%
  crop(.,hfbouds1)%>%
  projectRaster(.,crs='+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs')

hfDEM_matrix<-raster_to_matrix(hfDEM)
hfDEM_matrix %>%
  sphere_shade(texture ="imhof1")%>%
  add_water(detect_water(hfDEM_matrix),color="imhof1")%>%
  plot_3d(hfDEM_matrix,zscale = 1000,fov = 70,theta = 15,zoom = 0.60,phi = 30,windowsize = c(1000,800))




# writeRaster(hfDEM,"G:/Rdata/Anhui/DEM/hfDEM.tif")
# writeRaster(hfDEM,"G:/Rdata/Anhui/DEM/hfDEM1.tif")


hfDEM<- raster("G:\\Rdata\\Anhui\\DEM\\hfDEM.tif")
#check they have the same Coordinate Reference System(CRS)crs(listlandsat)
crs(hfbouds)
crs(hfDEM)
crs(listlandsat)



hflandsat_matrixr<- rayshader::raster_to_matrix(hflandsat$LC08_L2SP_121038_20201024_20201105_02_T1_SR_B4)
hflandsat_matrixg<- rayshader::raster_to_matrix(hflandsat$LC08_L2SP_121038_20201024_20201105_02_T1_SR_B3)
hflandsat_matrixb<- rayshader::raster_to_matrix(hflandsat$LC08_L2SP_121038_20201024_20201105_02_T1_SR_B2)
hfDEM_matrix <-rayshader::raster_to_matrix(hfDEM)


# rgb_array <- array(0,dim=c(nrow(hflandsat_matrixr),ncol(hflandsat_matrixr),3))
# rgb_array[,,1]<-hflandsat_matrixr/255
# rgb_array[,,2]<-hflandsat_matrixg/255
# rgb_array[,,3]<-hflandsat_matrixb/255
# 
# rgb_array <- aperm(rgb_array,c(2,1,3))
# plot_map(rgb_array)
# rgb_contrast <- scales::rescale(rgb_array,to=c(0,1))
# plot(rgb_contrast,stretch='lin')


hfDEM_matrix %>%
  sphere_shade(texture = "hefei") %>%
  plot_map()