裁剪 netcdf 文件并绘图

Question

我正在尝试按月为给定位置绘制混合层深度。

数据文件：http://www.ifremer.fr/cerweb/deboyer/mld/Surface_Mixed_Layer_Depth.php 它是页面上的最后一个文件，但现在任何文件都可以使用。

remove(list=ls())
library(raster)


mld <- brick("/Users/mld_DReqDTm02_c1m_reg2.0.nc", stopIfNotEqualSpaced = FALSE, varname = "mld")
print(mld)
extent(mld) <- extent(0, 360, -90, 90)
mld[mld > 1e4] <- NA
mld180 <- rotate(mld)
names(mld180) = month.abb

pprj <- "+proj=laea +lat_0=-90 +lon_0=180 +datum=WGS84 +ellps=WGS84 +no_defs +towgs84=0, 0, 0"

这是我裁剪数据的地方，但我不确定我是否正确执行。当我需要给它宽纬度和经度时，它可以工作，但我只想给它一个坐标，南纬 61 度到南纬 65 度，东经 140 度，我得到一个错误，它不会绘图。

g4 <- rgeos::gBuffer(SpatialPoints(cbind(0, 0), proj4string = CRS(pprj)), 
                 width = spDists(rbind(c(-140, -65), c(-140, -61)), longlat = TRUE, segments = T) * 1000, 
                 quadsegs = 180)



target <- projectExtent(mld180, pprj)
 Warning message:
 In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
 48 projected point(s) not finite



mld_trans <- crop(projectRaster(mld180, target, CRS(pprj)), g4)
mld_trans <- mask(mld_trans, g4)

boxplot(mld_trans,las = 1, xlab="Month", ylab="MLD (m)")

我做错了什么，如何将数据裁剪到这个位置，然后按月绘制深度？

Answer 1

下载文件，打开并修复它（按照您的示例）

url <- "http://www.ifremer.fr/cerweb/deboyer/data/mld_DReqDTm02_c1m_reg2.0.nc"
download.file(url, basename(url), mode="wb")  
library(raster)
mld <- brick("mld_DReqDTm02_c1m_reg2.0.nc", stopIfNotEqualSpaced = FALSE, varname = "mld")
extent(mld) <- extent(0, 360, -90, 90)
mld[mld > 1e4] <- NA
mld180 <- rotate(mld)
names(mld180) = month.abb

最简单的裁剪方法是创建一个范围。请注意，x 的最小值和最大值不能都是 -140。鉴于数据具有 2 度的空间分辨率，您会得到一列 -140 到 -138 或 -142 到 -140。我还调整了 y 限制（你会得到相同的 (-65, -61) --- 因为裁剪将值捕捉到网格单元格边界）

e <- extent(-140, -138, -64, -60)
x <- crop(mld180, e)

如果你想使用另一个 crs，你现在可以这样做

pprj <- "+proj=laea +lat_0=-90 +lon_0=180 +datum=WGS84"
mld_trans <- projectRaster(x, crs=pprj, method="ngb")