#快捷键: Ctrl+Shift+M, %>%
%>% 又叫管道符函数,该函数能将%>% 左侧运行的结果传递给%>%右侧的函数,默认作为右侧函数的第一个参数的实际值。如果想传递给右侧函数的指定参数,则把右侧函数的待传递的参数的实际值写为"."即可,"."就指代了%>% 左侧运行的结果。
1.安装并调用相关包
library(tidyverse)#处理数据
library(sf)#矢量数据的地理计算
library(raster)#栅格数据的地理计算
2.raster读取tif/nc文件并画图
# 用raster读取一个tif/nc 文件
raster("D:/data/merra2 DUST nc_tiff rasterdata/198001.tif") -> dust
dust
plot(dust)
plot结果
3.裁剪需要的研究区域
# 把研究区的数据裁剪出来
#先导入研究区矢量shp
read_sf("D:/文件备份lishouweidata/I QSF data整理/矢量边界/全国行政区划矢量数据_WGS84/市界_region.shp") -> prov
prov %>%
st_union() %>% #将省等合并
st_sf() -> cnpoly#转化为sf数据,中国的多边形数据
cnpoly %>%
st_simplify(dTolerance = 0.1) %>% #简化,加快运行
plot()
dust %>%
crop(cnpoly) %>%
mask(cnpoly) -> cndust
plot(cndust)
4.批量循环
dir.create(Asia_tif)#创建新文件夹目录
fs::dir_ls()可以查看文件夹中的文件名
正则表达式提取文件名进行循环
该处如何命名遇到困难,总报错~~~
str_extract(x, "\\d{6}") -> month %>%
paste0("D:/data/Asia_tif",month,".tif") -> destfile#生成文件名
这两句不对。
全部代码如下
library(tidyverse)#处理数据
library(sf)#矢量数据的地理计算
library(raster)#栅格数据的地理计算
# 用raster读取一个tif/nc 文件
raster("D:/data/merra2_dust original data/M2TMNXAER.5.12.4ːMERRA2_100.tavgM_2d_aer_Nx.198001.nc4") -> dust
dust
plot(dust)
# 把研究区的数据裁剪出来
#先导入研究区矢量shp
read_sf("D:/文件备份lishouweidata/I QSF data整理/矢量边界/全国行政区划矢量数据_WGS84/市界_region.shp") -> prov
prov %>%
st_union() %>% #将省等合并
st_sf() -> Asiapoly#转化为sf数据,中国的多边形数据
Asiapoly %>%
st_simplify(dTolerance = 0.1) %>% #简化,加快运行
plot()
dust %>%
crop(Asiapoly) %>% #裁剪
mask(Asiapoly) -> cndust#掩膜
plot(cndust)
#循环
dir.create("D:/data/Asia_tif")#创建新文件夹目录
#读取选中的文件夹中的列表:f in fs::dir_ls()
for (f in fs::dir_ls("D:/data/merra2_dust original data")){
f %>%
str_extract(x, "\\d{6}") -> month %>%
paste0("D:/data/Asia_tif",month,".tif") -> destfile#生成文件名
#为了防止电脑停电、死机等重复运行,加一个判断,如果文件已经存在则不需要再次运行
if(!file.exists(destfile)){
raster(f) %>%
crop(Asiapoly)%>%#裁剪
mask(Asiapoly)%>%#掩膜
writeRaster("")#保存
}
}
lapply(fs::dir_ls("D:/data/merra2 DUST nc_tiff rasterdata"), function(f){
f %>%
str_extract(x, "\\d{6}") -> month %>%
paste0("D:/data/Asia_tif",month,".tif") -> destfile#生成文件名
#为了防止电脑停电、死机等重复运行,加一个判断,如果文件已经存在则不需要再次运行
if(!file.exists(destfile)){
raster(f) %>%
crop(Asiapoly)%>%#裁剪
mask(Asiapoly)%>%#掩膜
writeRaster("")#保存
}
})
经李哥指导成功了!!(笔记本端运行的,机房断网了)
所以现在只需要找到标准的亚洲shp然后就能完成亚洲1980-2023年各月份dust分布的tiff文件了!
需要注意的地方 :raster(f, varname = "DUSMASS") %>% #需要指定图层
f %>%
str_extract(f,"\\d{6}") -> month
#%>% 表示把 前面的传递给后面函数的第一个参数,所以 f 写重复了
全部代码记录!!(20230826)
library(tidyverse)#处理数据
library(sf)#矢量数据的地理计算
library(raster)#栅格数据的地理计算
# 用raster读取一个tif/nc 文件
raster("E:/R 学习与数据处理/Merra2 官网粉尘原始数据/M2TMNXAER.5.12.4ːMERRA2_100.tavgM_2d_aer_Nx.198001.nc4",varname ="DUSMASS") -> dust
dust
plot(dust)
# 把研究区的数据裁剪出来
#先导入研究区矢量shp
read_sf("E:/I数据/矢量边界/全国行政区划矢量数据_WGS84/市界_region.shp") -> prov
prov %>%
st_union() %>% #将省等合并
st_sf() -> Asiapoly#转化为sf数据,中国的多边形数据
Asiapoly %>%
st_simplify(dTolerance = 0.1) %>% #简化,加快运行
plot()
dust %>%
crop(Asiapoly) %>% #裁剪
mask(Asiapoly) -> cndust#掩膜
plot(cndust)
#循环
dir.create("D:/data/Asia_tif/")#创建新文件夹目录
#读取选中的文件夹中的列表:f in fs::dir_ls()
fs::dir_ls("E:/R 学习与数据处理/Merra2 官网粉尘原始数据") -> flst
flst[1] %>%
str_extract("\\d{6}") -> month
paste0("D:/data/Asia_tif/", month, ".tif") -> destfile
raster(flst[1],varname ="DUSMASS") %>%
crop(Asiapoly) %>%
mask(Asiapoly) %>%
writeRaster(destfile,overwrite=TRUE)
#lapply循环
lapply(fs::dir_ls("E:/R 学习与数据处理/Merra2 官网粉尘原始数据"), function(f){
print(f)
f %>%
str_extract("\\d{6}") -> month
paste0("D:/data/Asia_tif/",month,".tif") -> destfile # 生成文件名
# 为了防止重复运行,加一个判断,如果文件已经存在则不需要再次运行
if(!file.exists(destfile)){
raster(f, varname = "DUSMASS") %>%
crop(Asiapoly) %>% # 裁剪
mask(Asiapoly) %>% # 掩膜
writeRaster(destfile) # 保存
}
})
#输出画图检查
raster("D:/data/Asia_tif/198002.tif") -> rst
plot(rst)
raster("D:/data/Asia_tif/198005.tif") -> rst
plot(rst)