绘制 NetCDF 的部分并覆盖形状文件-R

Posted 2023-03-23

【中文标题】绘制 NetCDF 的部分并覆盖形状文件-R

【英文标题】：Plot section of a NetCDF and overlay a shape file- R

【发布时间】：2018-07-05 03:28:35

【问题描述】：

我一直在绘制一个 NetCDF 文件并覆盖一个形状文件。然而，这是一个全球性的形象，我只想要其中的一部分，即只有墨西哥。 我使用This tutorial 来绘制所有内容并且效果很好。但是我应该在哪个部分裁剪地图。 这是我的代码：

require(utils)
require(colorRamps)
require(RNetCDF)
require(rasterVis)
require(rgdal)


Datos <- open.nc("./ERAII_viento_2014_300.nc")

u <- var.get.nc(Datos, "u")
v <- var.get.nc(Datos, "v")
lon <- var.get.nc(Datos, "longitude")
lat <-- var.get.nc(Datos, "latitude")

u[u == -32767] <- NA
v[v == -32767] <- NA

u9 <- raster(t(u[ , ,9])[ncol(u):1, ])
v9 <- raster(t(v[ , ,9])[ncol(v):1, ])

y <- brick(u9,v9)
projection(y) <- CRS("+init=epsg:4326")
extent(y) <- c(min(lon), max(lon), min(lat), max(lat))


slope <- sqrt(y[[1]]^2 + y[[2]]^2)
aspect <- atan2(y[[1]], y[[2]])


cntry <- readOGR(dsn="./shape_countries", layer="country")

we <- crop(y, extent(c(0, (180* res(y)[1]), min(lat), max(lat))))
ww <- crop(y, extent(c((180 * res(y)[1]), 360 * res(y)[1], min(lat), max(lat))))
extent(ww) <- c(-180 * res(y)[1], 0, min(lat), max(lat))

#extent(ww) #-180.50:0 y -90:90
#extent(we) #0:180.5 y -90:90
y2 <- merge(ww,we)

slope2 <- sqrt(y2[[1]]^2 + y2[[2]]^2)
vectorplot(y2 * 1.5, isField = "dXY", region = slope2, margin = FALSE, par.settings = rasterTheme(region = matlab.like(n = 10)),narrows = 10000, at = -15:50) + layer(sp.polygons(cntry))

我从http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/ 下载了风数据注意：我使用nco.exe 将变量u 和v 从“short”转换为“float”。 我不记得我从哪里得到 shapefile，但我认为它来自 ArcGIS。

【参考方案1】：

不幸的是，该教程使事情变得比需要的复杂很多。您需要做的就是：

library(raster)
library(ncdf4)    

u <- brick("./ERAII_viento_2014_300.nc", var="u")
v <- brick("./ERAII_viento_2014_300.nc", var="v")

将数据裁剪到墨西哥：

mex <- getData('GADM', country='MEX', level=0)
ux <- crop(u, mex)
uv <- crop(v, mex)

我没有你的文件，但我下载了一个：

f <- "_grib2netcdf.nc"
b <- brick(f, var="u10")

请注意，经度从 0 到 360（气候学家这样做），而不是标准的 -180 到 180。

extent(b)
#class       : Extent 
#xmin        : -0.375 
#xmax        : 359.625 
#ymin        : -90.375 
#ymax        : 90.375 

解决这个问题：

bb <- rotate(b)

bb
#class       : RasterBrick 
#dimensions  : 241, 480, 115680, 31  (nrow, ncol, ncell, nlayers)
#resolution  : 0.75, 0.75  (x, y)
#extent      : -179.625, 180.375, -90.375, 90.375  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#data source : in memory
#names       : X2016.01.01, X2016.01.02, X2016.01.03, X2016.01.04, X2016.01.05, X2016.01.06, X2016.01.07, X2016.01.08, X2016.01.09, X2016.01.10, X2016.01.11, X2016.01.12, X2016.01.13, X2016.01.14, X2016.01.15, ... 
#min values  :   -17.75753,   -30.02454,   -27.94626,   -18.51657,   -24.09466,   -22.60716,   -26.89947,   -24.13576,   -19.69050,   -20.82237,   -29.23778,   -20.65125,   -20.09774,   -26.92337,   -27.70344, ... 
#max values  :    21.46595,    22.34927,    25.07379,    24.71817,    20.92774,    21.60935,    24.50690,    23.31480,    24.46388,    25.14835,    26.64731,    21.74223,    20.15436,    29.76378,    27.18361, ... 
#Date/time   : 2016-01-01, 2016-01-31 (min, max)

现在你可以这样做了：

mex <- getData('GADM', country='MEX', level=0)
bbmex <- crop(bb, mex)

【讨论】：

谢谢，我确实得到了墨西哥的形状，但是当我使用“crop”时，我得到了这个错误： .local(x, y, ...) 中的错误：范围不重叠。我读到您可以使用“范围”功能来查看重叠是否存在，或者我做错了什么。