亚洲1980-2023年各月份dust分布的tiff(批量将nc/tif文件裁剪)

#快捷键： 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)