属性和位置分离的nc如何转tif

大多数nc文件可以直接使用terra包的 rast 函数读取，不过也有一些nc文件的属性信息和位置信息是分离的，直接读取的位置是错误的，就需要一些额外的处理。

实验数据获取

接下来以TAP大气污染物nc数据为例，介绍一下这种属性和位置分离的nc如何转为tif。

• TAP大气污染数据获取网址：http://tapdata.org.cn/?page_id=1162&item=pm25
• 从这里申请的月均NetCDF数据就是属性和位置分开的

TAP提供的nc数据如下所示，由两部分组成，一个是Grid_lonlat.nc，这个存储了位置信息，其他的nc里面是对应的PM2.5颗粒物属性信息。

解决思路

1. 读取两种nc文件，找到属性和位置信息之间的对应关系
2. 栅格转数据框，将属性和位置信息挂接起来
3. 将数据框转为矢量点
4. 矢量点栅格化
5. 输出转换好的tif文件

具体操作

数据读取和检查

读取两种nc文件，预览他们的缩略图和元数据，查看两个nc文件里面都有什么内容，找到他们之间的关联。

library(terra)
library(tidyverse)

PMlist = list.files(path = "./PM/China_PM25_1km_2023_01.nc~1/", pattern = ".nc$")
PMdir = paste0("./PM/China_PM25_1km_2023_01.nc~1/", PMlist)
PM25data = rast(PMdir[2])
PM25data
plot(PM25data)
pm25lonlat = rast(PMdir[7])
plot(pm25lonlat)
pm25lonlat

从上面可以发现，两个nc文件的行列一致，但是波段信息不一样，污染物数据属性所在的PM25data中只有一个波段，存储了PM2.5浓度信息，位置nc中存在两个波段，分别存储了Longitude经度和Latitude纬度信息。

位置和属性的关联

两个nc之间行列是一致的，因此我们可以使用像元位置xy将经纬度和PM2.5浓度属性联系起来，为保证挂接正确，我将XY坐标连起来作为唯一的标示。代码如下：

#栅格转数据框，添加像元坐标，为保证挂接正确，将XY坐标连起来作为唯一的标示
PM25_df = terra::as.data.frame(PM25data, xy= T) %>% 
  mutate(xy = str_c(x, y, sep = ""))
pm25lonlat_df = terra::as.data.frame(pm25lonlat, xy = T) %>% 
  mutate(xy = str_c(x, y, sep = ""))

#将两个nc合并成一个数据框
pm25nc_df = PM25_df %>% 
  dplyr::left_join(pm25lonlat_df, by = "xy") %>% 
  dplyr::select(c("Longitude", "Latitude", "PM25")) %>% 
  set_names(c("x", "y", "PM25"))

转栅格并输出

接下来把数据框转为了矢量点，然后根据矢量点的空间范围创建了一个基准栅格，然后利用这个基准栅格将矢量栅格化，得到属性、位置正确组合后的新栅格。这样就完成了两个nc之间信息的重组，输出后得到了属性和位置均正确的TIF数据。

#先生成一个矢量
pm25vect = vect(pm25nc_df, geom = c("x", "y"), crs = "EPSG:4326")

#生成一个标准栅格，然后将矢量栅格化
#根据矢量四至，经纬度nc行列号创建一个基准栅格
rast_template = rast(ext(pm25vect), 
                     nrows = ncol(pm25lonlat), ncols = nrow(pm25lonlat),
                     crs = "EPSG:4326")
rast_template
#将矢量栅格化
pm25rast = rasterize(pm25vect, rast_template, field = "PM25")
plot(pm25rast)
writeRaster(pm25rast, "./PM/PM25rast2302.tif")

自定义函数批处理

上面给出的代码是将一个nc进行了重组，有多套这样的nc数据时，需要使用批处理的方法来批量完成属性和空间位置的重组工作。

在R语言中，可以通过定义函数和循环的方式，来完成批处理。由于栅格的计算量比较大，对电脑内存消耗较大，因此没有使用并行计算。

pm25lonlat = rast(PMdir[7])
#创建一个自定义函数，输入nc文件名
fun_nctotif = function(inputnc){
  PM25data = rast(paste0("./PM/China_PM25_1km_2023_01.nc~1/", inputnc))
  PM25_df = terra::as.data.frame(PM25data, xy= T) %>% 
    mutate(xy = str_c(x, y, sep = ""))
  pm25lonlat_df = terra::as.data.frame(pm25lonlat, xy = T) %>% 
    mutate(xy = str_c(x, y, sep = ""))
  
  #将两个nc合并成一个数据框
  pm25nc_df = PM25_df %>% 
    dplyr::left_join(pm25lonlat_df, by = "xy") %>% 
    dplyr::select(c("Longitude", "Latitude", "PM25")) %>% 
    set_names(c("x", "y", "PM25"))
  
  #先生成一个矢量
  pm25vect = vect(pm25nc_df, geom = c("x", "y"), crs = "EPSG:4326")
  
  #生成一个标准栅格，然后将矢量栅格化
  #根据矢量四至，经纬度nc行列号创建一个标准栅格
  rast_template = rast(ext(pm25vect), 
                       nrows = ncol(pm25lonlat), ncols = nrow(pm25lonlat),
                       crs = "EPSG:4326")
  #将矢量栅格化
  pm25rast = rasterize(pm25vect, rast_template, field = "PM25")
  plot(pm25rast)
  writeRaster(pm25rast, paste0("./PM/", str_remove(inputnc, "\\.nc$"), ".tif"))
}

#通过循环完成批处理
for (i in 1:6) {
  fun_nctotif(PMlist[i])
}

答疑：多SHP多TIF情况下的提取栅格值到多个SHP

提问：有多个SHP，多个栅格的情况下，如何将多个栅格中的像元值，提取到每个SHP对应的点位？

回答：这种情况下可以使用嵌套循环，将SHP和TIF分别循环。代码如下， 详细讲解见阅读原文学习同名课程。

#读取多波段栅格
prelist = list.files(path = "./preYearsum2012-2022/", pattern = ".tif$")
predir = paste0("./preYearsum2012-2022/", prelist)
predir
pre = rast(predir)
plot(pre)
pre

#读取多个SHP
shplist = list.files(path = "./shp/", pattern = ".shp$")
shpdir = paste0("./shp/", shplist)
shpdir

#提取值到点
for (i in 1:length(shpdir)) {
  pt = vect(shpdir[i])
  for (j in 1:length(predir)) {
    pre = rast(predir[j])
    pre_pt = terra::extract(pre, pt, xy = T)
    write.csv(pre_pt, file = paste0("./csv/", str_remove(prelist[j], "\\.tif$"),
                                    str_remove(shplist[i], "\\.shp$"), ".csv"))
  }
}


paste0("./csv/", str_remove(prelist[1], "\\.tif$"), str_remove(shplist[1], "\\.shp$"))

str_remove(shplist[1], "\\.shp$")

参考文献

1. R语言栅格像元值提取，分区统计，提取采样点对应的像元值，栅格转数据框
2. nc文件terra包读取位置显示错误如何处理
3. NetCDF和HDF数据转TIFF，单波段和多波段栅格转换输出