作者：黄天元，



本帖子会简单介绍如何读入并处理栅格数据。首先，我们会用到一个矢量数据，数据来自：https://gadm.org/download_country_v3.html，用到的是澳洲的地图。读取方法如下：



# 获得数据的方法之一
# wget --no-check-certificate https://biogeo.ucdavis.edu/data/gadm3.6/gpkg/gadm36_AUS_gpkg.zip

library(pacman)
p_load(sf,raster,tidyverse)

# 查看有哪些图层
st_layers(
  "data/gadm36_AUS.gpkg"
)

# 读取特定图层
Ausoutline <- st_read("data/gadm36_AUS.gpkg", 
                      layer='gadm36_AUS_0')





R语言如何读取矢量数据 r语言读取tif_envi栅格TIF数据进行分割


可以对这个数据集进行观察：


# 观察
print(Ausoutline)

# 查看proj4坐标系
st_crs(Ausoutline)$proj4string


R语言如何读取矢量数据 r语言读取tif_envi栅格TIF数据进行分割_02


plot(Ausoutline)


R语言如何读取矢量数据 r语言读取tif_数据_03


然后，让我们读入栅格数据，这是一个全球平均气温数据，来源于：Historical climate data（tavg 5m）。


# 读入栅格数据
jan<-raster( "data/wc2.1_5m_tavg/wc2.1_5m_tavg_01.tif")
# have a look at the raster layer jan
jan
plot(jan)


R语言如何读取矢量数据 r语言读取tif_插值_04


R语言如何读取矢量数据 r语言读取tif_数据_05


我们可以改变它的坐标系，然后再次绘图：


# set the proj 4 to a new object
newproj<-"+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# get the jan raster and give it the new proj4
pr1 <- jan %>%
  projectRaster(., crs=newproj)
plot(pr1)


R语言如何读取矢量数据 r语言读取tif_R语言如何读取矢量数据_06


下面，让我们一次性读入所有栅格数据：


# 一次性读入所有栅格数据
dir("data/wc2.1_5m_tavg",full.names = T) %>% 
  stack() -> worldclimtemp
worldclimtemp


R语言如何读取矢量数据 r语言读取tif_envi栅格TIF数据进行分割_07


让我们为每一个图层命名，它其实是每个月的平均温度，因此可以用月份命名：


month <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", 
           "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

names(worldclimtemp) <- month


如果想取出一月份的数据，有两种方法：


worldclimtemp[[1]]
worldclimtemp$Jan


下面，让我们提取特定城市的数据，这是通过经纬度来提取的：


## 提取特定城市的数据
site <- c("Brisbane", "Melbourne", "Perth", "Sydney", "Broome", "Darwin", "Orange", 
          "Bunbury", "Cairns", "Adelaide", "Gold Coast", "Canberra", "Newcastle", 
          "Wollongong", "Logan City" )
lon <- c(153.03, 144.96, 115.86, 151.21, 122.23, 130.84, 149.10, 115.64, 145.77, 
         138.6, 153.43, 149.13, 151.78, 150.89, 153.12)
lat <- c(-27.47, -37.91, -31.95, -33.87, 17.96, -12.46, -33.28, -33.33, -16.92, 
         -34.93, -28, -35.28, -32.93, -34.42, -27.64)
#Put all of this inforamtion into one list 
samples <- data.frame(site, lon, lat, row.names="site")
samples
# Extract the data from the Rasterstack for all points 
AUcitytemp<- raster::extract(worldclimtemp, samples)


R语言如何读取矢量数据 r语言读取tif_数据_08


提取的变量是一个矩阵，可以稍微进行转化，把城市名称也加上。


AUcitytemp
Aucitytemp2 <- AUcitytemp %>% 
  as_tibble()%>% 
  add_column(Site = site, .before = "Jan")
Aucitytemp2


R语言如何读取矢量数据 r语言读取tif_R语言如何读取矢量数据_09


R语言如何读取矢量数据 r语言读取tif_插值_10


现在，我们要从世界气温图中取出澳洲的部分，实现如下（利用crop函数）：


Austemp <- Ausoutline %>%
  # now crop our temp data to the extent
  crop(worldclimtemp,.)

# plot the output
plot(Austemp)


R语言如何读取矢量数据 r语言读取tif_数据_11


这个图中还包括了非澳洲的部分，如果只想要澳洲的部分，可以这样操作：


exactAus <- Austemp %>%
  mask(Ausoutline, na.rm=TRUE)
plot(exactAus)


R语言如何读取矢量数据 r语言读取tif_envi栅格TIF数据进行分割_12


可以尝试对三月份的数据做一个温度分布直方图：


exactAusdf <- exactAus %>%
  as.data.frame()
# set up the basic histogram
gghist <- ggplot(exactAusdf, 
                 aes(x=Mar)) + 
  geom_histogram(color="black", 
                 fill="white")+
  labs(title="Ggplot2 histogram of Australian March temperatures", 
       x="Temperature", 
       y="Frequency")
# add a vertical line to the hisogram showing mean tempearture
gghist + geom_vline(aes(xintercept=mean(Mar, 
                                        na.rm=TRUE)),
                    color="blue", 
                    linetype="dashed", 
                    size=1)+
  theme(plot.title = element_text(hjust = 0.5))


R语言如何读取矢量数据 r语言读取tif_ci_13


更多可视化的选择可参考：https://andrewmaclachlan.github.io/CASA0005repo/rasters-descriptive-statistics-and-interpolation.html#histogram-with-ggplot

最后，我们会演示一下空间插值的过程。首先，我们来合并气温的时空分布数据：


samplestemp<-AUcitytemp%>%
  cbind(samples)
samplestemp


R语言如何读取矢量数据 r语言读取tif_数据_14


而后，让我们把它转换为空间信息数据框，这里需要声明经纬所在列，以及坐标系（这里定义为4326，也就是WGS 84）：


samplestemp<-samplestemp%>%
  st_as_sf(.,coords = c("lon", "lat"), 
           crs =4326, 
           agr = "constant")
samplestemp


R语言如何读取矢量数据 r语言读取tif_ci_15


这里我们可以对澳大利亚轮廓图和城市分布进行可视化：


plot(Ausoutline$geom)
plot(st_geometry(samplestemp), add=TRUE)


R语言如何读取矢量数据 r语言读取tif_ci_16


插值之前，让我们先统一坐标系，然后转换为SP对象：


samplestemp <- samplestemp %>%
  st_transform(3112)
Ausoutline <- Ausoutline %>%
  st_transform(3112)
samplestempSP <- samplestemp %>%
  as(., 'Spatial')
AusoutlineSP <- Ausoutline %>%
  as(., 'Spatial')


现在意味着，我们要用手头若干个城市的数值，来对整个澳大利亚的气温空间分布进行插值。我们需要创建一个要插值的空间范围：


emptygrd <- as.data.frame(spsample(AusoutlineSP, n=1000, type="regular", cellsize=200000))

names(emptygrd) <- c("X", "Y")

coordinates(emptygrd) <- c("X", "Y")

gridded(emptygrd) <- TRUE  # Create SpatialPixel object
fullgrid(emptygrd) <- TRUE  # Create SpatialGrid object

# Add the projection to the grid
proj4string(emptygrd) <- proj4string(samplestempSP)


然后进行插值：


p_load(gstat)
# Interpolate the grid cells using a power value of 2 
interpolate <- gstat::idw(Jan ~ 1, samplestempSP, newdata=emptygrd, idp=2.0)


这里用的IDW算法来插值，只使用1月份的数据。idp参数越大，则随着距离的增大，数值减少越大。如果idp为0，那么随着距离加大，依然不会有任何数值衰减。关于方法细节，可参考：How inverse distance weighted interpolation works。最后，我们看一下插值之后的结果：


# Convert output to raster object 
ras <- raster(interpolate)
# Clip the raster to Australia outline
rasmask <- mask(ras, Ausoutline)
# Plot the raster
plot(rasmask)


R语言如何读取矢量数据 r语言读取tif_数据_17


对空间插值感兴趣的读者，可以移步：Interpolation in R | Geodesic geometry

参考资料：

https://andrewmaclachlan.github.io/CASA0005repo/rasters-descriptive-statistics-and-interpolation.html