R语言栅格数据批量归一化 r语言读取栅格数据

作者：黄天元

之前按照Chapter 7 Advanced raster analysis | CASA0005 Geographic Information Systems and Science中给出的代码，尝试对栅格数据预处理。结果从数据下载到预处理都碰到了不少的问题。案例中用了Landsat8的数据，因此需要对Landsat数据有一定了解。这里直接给出相关链接：

刘大圣：Landsat介绍来啦

遥感小能人：Landsat8波段介绍

数据比较大，所在在本地处理的时候需要有一些技巧。下面，我们先读取数据（数据获取方法见Chapter 7 Advanced raster analysis | CASA0005 Geographic Information Systems and Science）：

library(pacman)
p_load(sp,raster,rgeos,rgdal,rasterVis,raster,fs,sf,tidyverse)

# read data
manchester_boundary <- st_read("data/manchester_boundary_download/Manchester_boundary.shp")
listlandsat <- dir("data/LC08_L1TP_203023_20190513_20190521_01_T1",
    pattern = "[B123456790].TIF",full.names = T) %>% stack()
#check they have the same Coordinate Reference System (CRS)
crs(manchester_boundary)
crs(listlandsat)

这里，我们把遥感数据读入到listlandsat中，而曼彻斯特的边界地图读到了manchester_boundary变量中。不过，我们发现第8个波段的分辨率与其他不一致，因此没有放倒stack函数中，我们需要单独对它进行重采样（resample），让所有图层保持一致。关于这方面的知识介绍，见

Resample-ArcGIS Pro | 文档pro.arcgis.com

Resample functiondesktop.arcgis.com

数据读取成功后，我们需要先把遥感影像限制在较小范围（曼彻斯特），否则每次处理都会牵涉大量数据。这一步非常重要，否则会非常耗时。

# get only Manchester
lsatmask <- listlandsat %>%
  # now crop our temp data to the extent
  crop(.,manchester_boundary)%>%
  mask(.,  manchester_boundary)

然后，我们用raster函数单独拉出第8个波段，并让它按照第1个波段的分辨率进行采样。采样方法为近邻采样法（nearest neighbour sampling，缩写为ngb），关于更多图像插值方法，见Rocky X：来聊聊图像插值算法（线性插值与非线性插值）。

# handle 8th band
b8list = dir("data/LC08_L1TP_203023_20190513_20190521_01_T1",
    pattern = "[B8].TIF",full.names = T) %>% 
  raster()

## ngb is a nearest neighbour sampling method
b8correct <- b8list%>%
  # now crop our temp data to the extent
  crop(.,manchester_boundary)%>%
  mask(.,  manchester_boundary)%>%
  resample(., lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B1, 
           method = "ngb")

插值后，让我们把数据叠加到原来的图层中。

lsatmask <- lsatmask %>%
  addLayer(., b8correct)

最后，我们要读出数据。我们会利用图层原始名称作为文件名，但是在最后加上“mask”，表明是从中挖出来的。我们会建立一个mask文件夹，来放置这些文件。

# add mask to the filenames within the raster stack
names(lsatmask) <- names(lsatmask)%>%
  str_c(., 
        "mask", 
        sep="_")

# I need to write mine out in another location
outputfilenames <-
  str_c("data/LC08_L1TP_203023_20190513_20190521_01_T1/mask/", names(lsatmask) ,sep="")

lsatmask %>%
  writeRaster(., outputfilenames, 
              bylayer=TRUE, 
              format='GTiff', 
              overwrite=TRUE)

最后，生成的文件就会导出到文件夹中。如下图所示：

本次分享就到这里结束，下次将会利用这些数据，继续进行更多的栅格数据高级操作，敬请期待。