如何使用 R 将区域从 shapefile 传输到 netcdf 文件？

编程问答 2022-05-02

问题描述

我有一个包含叶绿素数据的区域 netcdf 文件，我想添加一个新变量“区域”来映射有关海洋区域的信息。所以对于每个坐标点，都会有这个点属于哪个区域的信息。

因此，对于位于地中海地区的 netcdf 文件中的所有点，变量“region”将包含值 0，对于位于北海地区的所有点，变量“region”将包含值例如将包含值 1 等等...

我发现了一个包含基于生物地球化学过程（朗赫斯特省）的海洋区域的 shapefile。现在我想使用这个 shapefile 在 netcdf 文件中定义我的“区域”变量。有人知道我如何在 R 中做到这一点吗？

我认为我的问题与 DKRZ 的这篇帖子很接近，但我不想只提取/屏蔽一个区域，而是映射 shapefile 中定义的所有区域。

文件可以在这里找到：

https://drive.google.com/file/d/1kgPpHFapmuclDyUvw2TH_10i018GE9YH/view?usp=sharing

https://drive.google.com/file/d/1WLYEUs6XllZv6i0xjRif-N0syhR2lX01/view?usp=sharing

非常感谢！

编辑
我发现这篇文章帮助我解决了我的问题： https://www.r-bloggers.com/2014/05/converting-shapefiles-to-rasters-in-r/)

解决方法

这是我可以根据 Map Lat/Lon Points to a Shape File in R 组合在一起的内容。主要工作由 over() 完成，需要几分钟时间。

library(ncdf4)
nc_data = nc_open('/mnt/d/Downloads/CCI_ALL-v5.0-MONTHLY.nc',T)    # write
lon = ncvar_get(nc_data,"lon")
lat = ncvar_get(nc_data,"lat")
library(rgdal)
shapefile = readOGR('/mnt/d/Downloads/Longhurst_world_v4_2010.shp')
grid = expand.grid(lon=lon,lat=lat)
coordinates(grid) = ~lon+lat
grid@proj4string = shapefile@proj4string
# for each grid point,find the province it belongs to (if available):
prov = over(grid,shapefile['ProvCode'])
# define dimensions for the new variable
londim = ncdim_def("lon",ncatt_get(nc_data,"lon")$units,lon)
latdim = ncdim_def("lat","lat")$units,lat)
# define the new variable
region = ncvar_def("region","",list(londim,latdim),prec="integer")
# add the new variable to the netCDF file
nc_data = ncvar_add(nc_data,region)
# write the variable's values to the file
ncvar_put(nc_data,region,prov$ProvCode)
nc_close(nc_data)

rgdal 版本等：

Loading required package: sp
rgdal: version: 1.4-8,(SVN revision 845)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 3.0.4,released 2020/01/28
 Path to GDAL shared files:
 GDAL binary built with GEOS: TRUE
 Loaded PROJ.4 runtime: Rel. 6.3.1,February 10th,2020,[PJ_VERSION: 630]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.3-2
OGR data source with driver: ESRI Shapefile
Source: "/mnt/d/Downloads/Longhurst_world_v4_2010.shp",layer: "Longhurst_world_v4_2010"
with 54 features
It has 2 fields
mappingmappingnetcdfoceanrrshapefile