envi栅格TIF数据进行分割_R语言空间数据分析（五）：栅格数据处理

作者：黄天元，复旦大学博士在读，热爱数据科学与开源工具（R），致力于利用数据科学迅速积累行业经验优势和科学知识发现，涉猎内容包括但不限于信息计量、机器学习、数据可视化、应用统计建模、知识图谱等，著有《R语言高效数据处理指南》（《R语言数据高效处理指南》(黄天元)【摘要书评试读】- 京东图书）。知乎专栏：R语言数据挖掘。邮箱：huang.tian-yuan@qq.com.欢迎合作交流。

HopeR：R语言空间数据分析（零）：总目录zhuanlan.zhihu.com

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

# 获得数据的方法之一
# wget --no-check-certificate https://biogeo.ucdavis.edu/data/gadm3.6/gpkg/gadm36_AUS_gpkg.ziplibrary(pacman)
p_load(sf,raster,tidyverse)# 查看有哪些图层
st_layers("data/gadm36_AUS.gpkg"
)# 读取特定图层
Ausoutline <- st_read("data/gadm36_AUS.gpkg", layer='gadm36_AUS_0')

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

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

plot(Ausoutline)

然后，让我们读入栅格数据，这是一个全球平均气温数据，来源于：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)

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

# 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)

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

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

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

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)

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

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

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

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

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

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

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

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))

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

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

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

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

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

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

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

插值之前，让我们先统一坐标系，然后转换为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)

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

参考资料：

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

envi栅格TIF数据进行分割_R语言空间数据分析（五）：栅格数据处理相关推荐

r语言worldclim数据_R语言空间数据分析（五）：栅格数据处理
作者:黄天元,复旦大学博士在读,热爱数据科学与开源工具(R),致力于利用数据科学迅速积累行业经验优势和科学知识发现,涉猎内容包括但不限于信息计量.机器学习.数据可视化.应用统计建模.知识图谱等,著有& ...
moran指数 r语言_R语言空间数据分析（七）：空间自相关
作者:黄天元,复旦大学博士在读,热爱数据科学与开源工具(R),致力于利用数据科学迅速积累行业经验优势和科学知识发现,涉猎内容包括但不限于信息计量.机器学习.数据可视化.应用统计建模.知识图谱等,著有& ...
envi栅格TIF数据进行分割_栅格数据镶嵌
试了很多坑,最后发现arcgis用来栅格数据的镶嵌,效果是不言而喻,接下来就简单介绍流程,希望对没有跳过坑的小伙伴有用(以下均用英文arcgis版本操作)!!! 1 首先打开arggis中ArcToo ...
envi栅格TIF数据进行分割_常用水文气象数据读取及其可视化（二进制、HDF5、NetCDF）以GLDAS、MODIS、GSMaP为例...
" 地学.水文.气象领域的自然科学数据通常以netcdf.hdf.二进制等方式存储,比如温度.降水.蒸发数据等:学会这些数据格式的读取和可视化是进行地学统计分析计算的关键,python提供了 ...
语言转4字节数据整型_R语言与RGui平台_数据类型_向量_4
计算机语言RGui平台上的R语言__数据类型_向量_4 R语言进阶4_数据类型_向量1 咱们而可以从自然语言(汉语.英.法.俄.德.日.拉丁.伊斯兰.等等)的基本特征来看--词汇.句子.段落.结构.文 ...
合并相同数据的行_R语言笔记（六）：数据框重塑(reshape2)
数据处理主要内容包括: 1. 特殊值处理 1.1 缺失值 1.2 离群值 1.3 日期 2. 数据转换(base vs. dplyr) 2.1 筛选(subset vs. filter/select/ ...
python 一组数据正态分布散点图_R语言入门之散点图
散点图 1. 简单散点图在R中有很多方式去绘制散点图,其中最基本的就是是用plot(x, y)函数,往期内容已经进行过详细讲解,这里就不赘述了,下面直接看实例图. # 简单散点图 attach(mt ...
python 分析两组数据的差异_R语言limma包差异基因分析（两组或两组以上）
使用limma包进行差异基因分析时,做最多的是两分类的,例如control组和disease组,但也会碰到按照序列进行的分组.这时,如果逐一使用两两比较求差异基因则略显复杂.其实开发limma包的大神 ...
r语言数据变量分段_R语言：统计多个数据框中分类变量各值的频数
导读查看文件.获取ID 读取数据表读取数据表一.查看文件.获取ID 关键参数: list.files(pattern="条件") # 根据条件获取文件名 strsplit(向 ...

envi栅格TIF数据进行分割_R语言空间数据分析（五）：栅格数据处理

envi栅格TIF数据进行分割_R语言空间数据分析（五）：栅格数据处理相关推荐

最新文章

热门文章