R语言：栅格数据处理

最新推荐文章于 2024-10-16 10:13:45 发布

树谷-胡老师

最新推荐文章于 2024-10-16 10:13:45 发布

阅读量2.1k

点赞数 1

分类专栏：教程文章标签： java dreamweaver html

原文链接：https://zhuanlan.zhihu.com/p/280968987

版权

教程专栏收录该内容

147 篇文章 203 订阅

订阅专栏

本帖子会简单介绍如何读入并处理栅格数据。首先，我们会用到一个矢量数据，数据来自：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(  "H:/shp/gadm36_KEN_gpkg/gadm36_KEN.gpkg")Kenoutline <- st_read("H:/shp/gadm36_KEN_gpkg/gadm36_KEN.gpkg",                       layer='gadm36_KEN_2')#打印这个数据print(Kenoutline)#查看其坐标系st_crs(Kenoutline)# Coordinate Reference System:#   User input: WGS 84 # wkt:#   GEOGCRS["WGS 84",#           DATUM["World Geodetic System 1984",#                 ELLIPSOID["WGS 84",6378137,298.257223563,#                           LENGTHUNIT["metre",1]]],#           PRIMEM["Greenwich",0,#                  ANGLEUNIT["degree",0.0174532925199433]],#           CS[ellipsoidal,2],#           AXIS["geodetic latitude (Lat)",north,#                ORDER[1],#                ANGLEUNIT["degree",0.0174532925199433]],#           AXIS["geodetic longitude (Lon)",east,#                ORDER[2],#                ANGLEUNIT["degree",0.0174532925199433]],#           USAGE[#             SCOPE["unknown"],#             AREA["World"],#             BBOX[-90,-180,90,180]],#           ID["EPSG",4326]]
plot(Kenoutline)

读取栅格数据

tempreture<-raster("H:/shp/wc2.1_10m_tavg/wc2.1_10m_tavg_01.tif")tempreture# class      : RasterLayer # dimensions : 1080, 2160, 2332800  (nrow, ncol, ncell)# resolution : 0.1666667, 0.1666667  (x, y)# extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)# crs        : +proj=longlat +datum=WGS84 +no_defs # source     : H:/shp/wc2.1_10m_tavg/wc2.1_10m_tavg_01.tif # names      : wc2.1_10m_tavg_01 # values     : -45.884, 34.0095  (min, max)st_crs(tempreture)# Coordinate Reference System:#   User input: +proj=longlat +datum=WGS84 +no_defs # wkt:#   GEOGCRS["unknown",#           DATUM["World Geodetic System 1984",#                 ELLIPSOID["WGS 84",6378137,298.257223563,#                           LENGTHUNIT["metre",1]],#                 ID["EPSG",6326]],#           PRIMEM["Greenwich",0,#                  ANGLEUNIT["degree",0.0174532925199433],#                  ID["EPSG",8901]],#           CS[ellipsoidal,2],#           AXIS["longitude",east,#                ORDER[1],#                ANGLEUNIT["degree",0.0174532925199433,#                          ID["EPSG",9122]]],#           AXIS["latitude",north,#                ORDER[2],#                ANGLEUNIT["degree",0.0174532925199433,#                          ID["EPSG",9122]]]]plot(tempreture)

改变栅格投影坐标系

newproj<-"+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"tempreture1<-tempreture%>%projectRaster(crs = newproj)plot(tempreture1)

％>％来自dplyr包的管道函数，其作用是将前一步的结果直接传参给下一步的函数，从而省略了中间的赋值步骤，可以大量减少内存中的对象，节省内存。

一次读入所有栅格数据

#一次读入所有栅格数据dir("G:/shp/wc2.1_10m_tavg",full.names = TRUE)%>%  stack()->worldclimateworldclimate# class      : RasterStack # dimensions : 1080, 2160, 2332800, 12  (nrow, ncol, ncell, nlayers)# resolution : 0.1666667, 0.1666667  (x, y)# extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)# crs        : +proj=longlat +datum=WGS84 +no_defs # names      : wc2.1_10m_tavg_01, wc2.1_10m_tavg_02, wc2.1_10m_tavg_03, wc2.1_10m_tavg_04, wc2.1_10m_tavg_05, wc2.1_10m_tavg_06, wc2.1_10m_tavg_07, wc2.1_10m_tavg_08, wc2.1_10m_tavg_09, wc2.1_10m_tavg_10, wc2.1_10m_tavg_11, wc2.1_10m_tavg_12 # min values :         -45.88400,         -44.80000,         -57.92575,         -64.19250,         -64.81150,         -64.35825,         -68.46075,         -66.52250,         -64.56325,         -55.90000,         -43.43475,         -45.32700 # max values :          34.00950,          32.82425,          32.90950,          34.19375,          36.25325,          38.35550,          39.54950,          38.43275,          35.79000,          32.65125,          32.78800,          32.82525
month<-c("Jan","Feb","Mar","Apr","May","Jun",         "Jul","Aug","Sep","Oct","Nov","Dec")names(worldclimate)<-monthworldclimateworldclimate$[[1]]worldclimate$Jan# class      : RasterStack # dimensions : 1080, 2160, 2332800, 12  (nrow, ncol, ncell, nlayers)# resolution : 0.1666667, 0.1666667  (x, y)# extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)# crs        : +proj=longlat +datum=WGS84 +no_defs # names      :       Jan,       Feb,       Mar,       Apr,       May,       Jun,       Jul,       Aug,       Sep,       Oct,       Nov,       Dec # min values : -45.88400, -44.80000, -57.92575, -64.19250, -64.81150, -64.35825, -68.46075, -66.52250, -64.56325, -55.90000, -43.43475, -45.32700 # max values :  34.00950,  32.82425,  32.90950,  34.19375,  36.25325,  38.35550,  39.54950,  38.43275,  35.79000,  32.65125,  32.78800,  32.82525 
#提取特定城市的温度数据Kenoutline$NAME_1site<-c("Mombasa", "Nakuru", "Meru","Marsabit","Kakuma")lon<-c(39.655,36.8172,37.64561,38.01902,34.85441)lat<-c(-4.0412,-1.28962,0.05186,2.33916,3.71634)#put all of this information into one listsamples<-data.frame(site,lon,lat,row.names = "site")samples# Exctrct the data from the Rasterstack for all pointKenCityTemp<-raster::extract(worldclimate,samples)class(KenCityTemp)KenCityTemp# Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec# [1,] 27.19083 27.56361 28.08833 27.31250 25.74750 24.70750 23.78306 23.92917 24.62333 25.64444 26.59528 27.02778# [2,] 19.37425 19.99025 20.23750 19.65600 18.58350 17.18000 16.34900 16.65275 17.83250 19.15825 18.91300 18.67950# [3,] 17.17025 17.95600 18.43675 17.91850 17.32725 16.31475 15.75400 16.03425 17.35950 18.02175 17.15275 16.67075# [4,] 24.05100 24.76150 24.79825 23.70475 23.18525 22.18050 21.42425 21.66375 22.63575 23.05825 22.76675 23.13600# [5,] 27.90425 28.59675 28.63800 27.48875 27.30100 26.96225 26.41325 26.43175 27.30975 27.54900 27.32775 27.19650KenCityTemp2<- KenCityTemp %>%   as_tibble()%>%  add_column(Site=site,.before="Jan")KenCityTemp2# Site       Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec# <chr>    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>#   1 Mombasa   27.2  27.6  28.1  27.3  25.7  24.7  23.8  23.9  24.6  25.6  26.6  27.0# 2 Nakuru    19.4  20.0  20.2  19.7  18.6  17.2  16.3  16.7  17.8  19.2  18.9  18.7# 3 Meru      17.2  18.0  18.4  17.9  17.3  16.3  15.8  16.0  17.4  18.0  17.2  16.7# 4 Marsabit  24.1  24.8  24.8  23.7  23.2  22.2  21.4  21.7  22.6  23.1  22.8  23.1# 5 Kakuma    27.9  28.6  28.6  27.5  27.3  27.0  26.4  26.4  27.3  27.5  27.3  27.2

从世界气温图中裁剪出肯尼亚的部分，使用crop函数

Kentemp<-Kenoutline %>%  #now crop our temp data to the extent  crop(worldclimate,.)Kentemp# class      : RasterBrick # dimensions : 58, 49, 2842, 12  (nrow, ncol, ncell, nlayers)# resolution : 0.1666667, 0.1666667  (x, y)# extent     : 33.83333, 42, -4.666667, 5  (xmin, xmax, ymin, ymax)# crs        : +proj=longlat +datum=WGS84 +no_defs # source     : memory# names      :      Jan,      Feb,      Mar,      Apr,      May,      Jun,      Jul,      Aug,      Sep,      Oct,      Nov,      Dec # min values :  7.00550,  7.59825,  7.86850,  7.01375,  6.01350,  5.03775,  4.30050,  4.78400,  5.54800,  6.64600,  6.82000,  6.88075 # max values : 30.18350, 30.95450, 31.74950, 30.16625, 29.74000, 29.21300, 28.52475, 28.69950, 29.58900, 30.09900, 29.56600, 29.33750 #plot the outputplot(Kentemp)

裁剪出真正属于肯尼亚的部分

exactKen<-Kentemp %>%  mask(Kenoutline,na.rm=TRUE)plot(exactKen)

或者可以这样(一步到位)：

Kentemp<-Kenoutline %>%  #now crop our temp data to the extent  crop(worldclimate,.) %>%  mask(.,Kenoutline,na.rm = TRUE)Kentemp#plot the outputplot(Kentemp)```对肯尼亚八月份的数据做一个温度直方图
```r#可以尝试对8月份的数据做一个温度分布直方图Kentempdf<-Kentemp %>%  as.data.frame()#set up the basic histogramgghist<-ggplot(Kentempdf,               aes(x=Aug))+  geom_histogram(color="black",                 fill="white")+  labs(title="Ggplot2 histogram of Kenya August Temperatures",       x="Temperatures",       y="Frequency")#add a vertical line to the histogram showing mean temperaturegghist + geom_vline(aes(xintercept=mean(Aug,                                        na.rm=TRUE)),                    color="blue",                    linetype="dashed",                    size=1)+  theme(plot.title = element_text(hjust = 0.5))

空间插值之前，先合并气温的时空分布数据

#最后，我们会演示一下空间插值的过程。首先，我们来合并气温的时空分布数据：samplestemp<-KenCityTemp%>%  cbind(samples)samplestemp# Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec# Mombasa  27.19083 27.56361 28.08833 27.31250 25.74750 24.70750 23.78306 23.92917 24.62333 25.64444 26.59528 27.02778# Nakuru   19.37425 19.99025 20.23750 19.65600 18.58350 17.18000 16.34900 16.65275 17.83250 19.15825 18.91300 18.67950# Meru     17.17025 17.95600 18.43675 17.91850 17.32725 16.31475 15.75400 16.03425 17.35950 18.02175 17.15275 16.67075# Marsabit 24.05100 24.76150 24.79825 23.70475 23.18525 22.18050 21.42425 21.66375 22.63575 23.05825 22.76675 23.13600# Kakuma   27.90425 28.59675 28.63800 27.48875 27.30100 26.96225 26.41325 26.43175 27.30975 27.54900 27.32775 27.19650# lon      lat# Mombasa  39.65500 -4.04120# Nakuru   36.81720 -1.28962# Meru     37.64561  0.05186# Marsabit 38.01902  2.33916# Kakuma   34.85441  3.71634

然后把它转化为空间信息数据框，这里需要声明经纬所在列，以及坐标系定义为4326（WGS 84）

#然后把它转化为空间信息数据框，这里需要声明经纬所在列，以及坐标系定义为4326samplestemp<-samplestemp%>%  st_as_sf(.,coords=c("lon","lat"),           crs=4326,           agr="constant")#samplestempplot(Kenoutline$geometry)plot(st_geometry(samplestemp),add=TRUE)

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

samplestemp# Simple feature collection with 5 features and 12 fields# Attribute-geometry relationship: 12 constant, 0 aggregate, 0 identity# geometry type:  POINT# dimension:      XY# bbox:           xmin: 34.85441 ymin: -4.0412 xmax: 39.655 ymax: 3.71634# geographic CRS: WGS 84# Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec# Mombasa  27.19083 27.56361 28.08833 27.31250 25.74750 24.70750 23.78306 23.92917 24.62333 25.64444 26.59528 27.02778# Nakuru   19.37425 19.99025 20.23750 19.65600 18.58350 17.18000 16.34900 16.65275 17.83250 19.15825 18.91300 18.67950# Meru     17.17025 17.95600 18.43675 17.91850 17.32725 16.31475 15.75400 16.03425 17.35950 18.02175 17.15275 16.67075# Marsabit 24.05100 24.76150 24.79825 23.70475 23.18525 22.18050 21.42425 21.66375 22.63575 23.05825 22.76675 23.13600# Kakuma   27.90425 28.59675 28.63800 27.48875 27.30100 26.96225 26.41325 26.43175 27.30975 27.54900 27.32775 27.19650# geometry# Mombasa    POINT (39.655 -4.0412)# Nakuru   POINT (36.8172 -1.28962)# Meru     POINT (37.64561 0.05186)# Marsabit POINT (38.01902 2.33916)# Kakuma   POINT (34.85441 3.71634)samplestemp<-samplestemp%>%  st_transform(21097)samplestempKenoutline<-Kenoutline%>%  st_transform(21097)
samplestempsp<-samplestemp%>%  as(.,'Spatial')Kenoutlinesp<-Kenoutline%>%  as(.,'Spatial')#现在意味着我们要用手头若干个城市的数值，来对澳大利亚的气温空间分布进行插值，#我们需要创建一个要插值的空间范围emptygrd<-as.data.frame(spsample(Kenoutlinesp,n=1000,type = "regular",cellsize=60000))names(emptygrd)<-c("X", "Y")plot(emptygrd)coordinates(emptygrd)<-c("X", "Y")gridded(emptygrd)<-TRUE #create SpatialPixel objectfullgrid(emptygrd)<-TRUE #create SpatialGrid object
#add the projection to the gridproj4string(emptygrd)<-proj4string(samplestempsp)
#然后进行插值：p_load(gstat)#Interpolate the grid cells using a power value of 2interpolate<-gstat::idw(Jan~1, samplestempsp, newdata=emptygrd, idp=2.0)#这里使用IDM算法来插值，只使用1月份的数据，idp参数越大，则随着距离的增大，#数值减少越大。如果idp为0，那么随着距离加大，依然不会有任何数值衰减。#convert output to raster objectras<-raster(interpolate)%>%  mask(.,Kenoutline)plot(ras)