r 读取 NetCDF 并导出为 shapefile

Question

我需要使用 R 读取 NetCDF 文件并将每个时间步导出为平滑的多边形形状文件。 我有两个问题：平滑光栅和从 NC 文件正确投影导出到 shapefile。 输出是一个规则的网格，没有投影。

这是一个示例代码：

>NCFileName = MyncFile.nc
NCFile = open.ncdf(NCFileName) 
NCFile 
[1] "file CF_OUTPUT.nc has 6 dimensions:"
[1] "time   Size: 61"
[1] "height   Size: 8"
[1] "lat   Size: 185"
[1] "lon   Size: 64"
[1] "Time   Size: 61"
[1] "DateStrLen   Size: 19"
[1] "------------------------"
[1] "file CF_OUTPUT.nc has 20 variables:"
[1] "float temp[lon,lat,height,time]  Longname:Temperature Missval:1e+30"
[1] "float relh[lon,lat,height,time]  Longname:Relative Humidity Missval:1e+30"
[1] "float airm[lon,lat,height,time]  Longname:Air density Missval:1e+30"
[1] "float z[lon,lat,height,time]  Longname:Layer top altitude Missval:1e+30"
[1] "float ZH[lon,lat,height,time]  Longname:Layer top altitude Missval:1e+30"
[1] "float hlay[lon,lat,height,time]  Longname:Layer top altitude Missval:1e+30"
[1] "float PM10ant[lon,lat,height,time]  Longname:PM10ant Concentration Missval:1e+30"
[1] "float PM10bio[lon,lat,height,time]  Longname:PM10bio Concentration Missval:1e+30"
[1] "float PM10[lon,lat,height,time]  Longname:PM10 Concentration Missval:1e+30"
[1] "float PM25ant[lon,lat,height,time]  Longname:PM25ant Concentration Missval:1e+30"
[1] "float PM25bio[lon,lat,height,time]  Longname:PM25bio Concentration Missval:1e+30"
[1] "float PM25[lon,lat,height,time]  Longname:PM25 Concentration Missval:1e+30"
[1] "float C2H4[lon,lat,height,time]  Longname:C2H4 Concentration Missval:1e+30"
[1] "float CO[lon,lat,height,time]  Longname:CO Concentration Missval:1e+30"
[1] "float SO2[lon,lat,height,time]  Longname:SO2 Concentration Missval:1e+30"
[1] "float NO[lon,lat,height,time]  Longname:NO Concentration Missval:1e+30"
[1] "float NO2[lon,lat,height,time]  Longname:NO2 Concentration Missval:1e+30"
[1] "float O3[lon,lat,height,time]  Longname:O3 Concentration Missval:1e+30"
[1] "char Times[DateStrLen,Time]  Longname:Times Missval:NA"
[1] "float HGT[lon,lat,time]  Longname:Topography Missval:1e+30"

nc.a=get.var.ncdf(NCFile , varid = 'NO2', start=c(1,1,1,1), count=c(-1,-1,1,1))
Pol <- rasterToPolygons(raster(nc.a),dissolve = TRUE)
Pol
class       : SpatialPolygonsDataFrame 
features    : 11829 
extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
variables   : 1
names       :             layer 
min values  : 0.219758316874504 
max values  :  0.84041428565979 
writeOGR(Pol, dsn = getwd(), layer = 'testPol', driver = 'ESRI Shapefile', overwrite_layer = TRUE)

然而，我得到的是未投影的网格多边形。

更新： 按照@kakk11 和@RobertH 的回答，我能够解决部分问题。我仍然得到一个网格状的多边形，没有平滑。这是我到目前为止所做的： 正如@RobertH 建议的那样，我无法将变量直接提取到栅格中。所以我使用了“get.var.ncdf”，然后使用了“raster”：

NCFileName = 'MyncFile.nc'
NCFile = open.ncdf(NCFileName)
nc.a = get.var.ncdf(NCFile, varid = 'NO2', start=c(1,1,1,13), count=c(-1,-1,1,1))
nc.a = raster(nc.a)
# put in correct extent:
lat  = NCFile$dim$lat$vals
lon  = NCFile$dim$lon$vals
ExtentLat = range(lat)
ExtentLon = range(lon)
rm(lat,lon)

nc.a = flip(t(nc.a), direction='y')

# Give it lat/lon coords 
extent(nc.a) = c(ExtentLon,ExtentLat)

然后'cut'命令返回向量，所以我使用'ratser:reclassify'：

cuts = c(0,5,15,30,50)
classes <- cbind(cuts[1:length(cuts)-1],cuts[2:length(cuts)],cuts[2:length(cuts)])
nc.class <- reclassify(nc.a, classes)

然后我使用 'rasterToPolygons' 和 'dissolve=TRUE' 创建多边形：

pol <- rasterToPolygons(nc.class, dissolve = TRUE)
# set UTM projection:
WGS84_Projection = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(pol) <- CRS(WGS84_Projection)
writeOGR(pol, dsn = getwd(), layer = 'file' , driver = 'ESRI Shapefile', overwrite_layer = TRUE)

不过，所有这些都会创建多边形 shapefile，其中多边形不平滑，这是主要挑战。

在这方面可以使用一些帮助。

喜欢

Answer 1

您首先需要正确创建一个 RasterLayer，如下所示：

 r <- raster('MyncFile.nc', var='NO2')
 # or, to get all time steps at once
 # brick('MyncFile.nc', var='NO2')

然后您可以使用重新分类或切割来概括（分类）这些值。例如

 cuts <- seq(0.2, 0.9, 0.1)
 rc <- cut(r, cuts)

制作多边形并保存到 shapefile

 pol <- rasterToPolygons(rc, dissolve = TRUE)
 shapefile(pol, 'file.shp')