裁剪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, pprj)