r - 如何设置raster::plot的显示范围？

Question

摘要

如何将R raster::plot的显示限制为Raster *对象的边界？为什么我问：

我是R初学者，必须


将未投影的（经纬度）全局空间数据从ASCII转换为netCDF（已解决）
“重新定义”它：即，将数据限制在一个区域中，对其进行投影，然后将其缩小（减小网格单元的大小）（大部分已解决）


不幸的是，在科学工作中，主要是QA数据转换，例如通过视觉检查进行的转换。上面的步骤1看起来不错。但是尽管第2步现在看起来要好得多，但我确实只想绘制重新定义的extents的边界（或RasterLayer）。我在公共git存储库（here）中有由bash脚本驱动的R代码，该库执行转换步骤，并绘制转换后的输出apparently correctly。但是，我用raster::plot绘制重新网格的输出的尝试不太正确（请参见输出here的前三页）。

怎么修？

细节

背景

我需要获取一些数据（全球海洋排放清单（EI）），并将其与其他数据（主要是其他EI）结合起来，然后输入到模型中。该模型想要消耗netCDF，并且希望该netCDF在比连续美国（CONUS）稍大的空间域上使用。该图像的边界是域的边界。 （请注意，该域的一部分，但只有一部分是海洋的。）该模型还希望netCDF以某种方式投影（LCC的分辨率为12 km）。我将其视为2个独立的问题：


将全局EI从其本机ASCII格式转换为netCDF
将netCDF从全局/未投影“缩减”到缩小的（更精细的分辨率）投影子域。


我正在尝试使用我在this public git repository处存档的代码来解决这些问题。如果克隆该git repo，然后按照“第一个示例” in the README所述配置/运行bash驱动程序GEIA_to_netCDF.sh（并假定已正确配置R，尤其是使用packages = raster，ncdf4）它将输出netCDF并绘制（使用fields::image.plot）输出到PDF的图，该图应如下所示：。输出数据的分配似乎正确，因此问题1似乎已解决。

问题

解决问题2（又称“第二个示例” in the README）似乎需要R package = raster。我的bash驱动程序regrid_GEIA_netCDF.sh调用regrid.global.to.AQMEII.r（均为in repository），该命令运行无误，并显示其输出的PDF。 GEIA_N2O_oceanic_regrid.pdf包含4页，对应于regrid.global.to.AQMEII.r末尾的4个代码块。这里只有前3页有用（第4页尝试使用fields::image.plot并有更大的问题）。

Page = 1/4的结果来自简单的raster::plot重新生成的输出，加上来自wrld_simpl的CONUS地图的投影版本：

plot(out.raster)
plot(map.us.proj, add=TRUE)

不幸的是，输出看起来是全局的，或几乎是全局的：但是将输出重新定义到的期望域要小得多：。因此，我有3个问题（1个主要问题，2个后续问题）：


默认情况下，由 raster::plot显示的图像仅显示其参数中指定的 RasterLayer范围吗？
如果是这样，我的后悔何在？ （下面有更多内容）
如果不是（即默认情况下 raster::plot显示的内容多于其 RasterLayer的范围），如何使 raster::plot仅显示该 RasterLayer的范围？


regrid.global.to.AQMEII.r（ here）中进行重新网格化的代码似乎获得了正确的输出：

out.crs <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000'

请注意，CRS似乎与给定 here的输出域的定义匹配。 （在链接页面上，滚动到地图。）

out.raster <-
  projectRaster(
from=in.raster, to=template.raster, method='bilinear', crs=out.crs,
overwrite=TRUE, progress='window', format='CDF',
# args from writeRaster
NAflag=-999.0,  # match emi_n2o:missing_value,_FillValue (TODO: copy)
varname=data.var.name, 
varunit='ton N2O-N/yr',
longname='N2O emissions',
xname='COL',
yname='ROW',
filename=out.fp)
out.raster
# made with CRS
#> class       : RasterLayer 
#> dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
#> resolution  : 53369.55, 56883.69  (x, y)
#> extent      : -14802449, 9694173, -6258782, 10749443  (xmin, xmax, ymin, ymax)
# ??? why still proj=longlat ???
#> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#> data source : /home/rtd/code/R/GEIA_to_netCDF/GEIA_N2O_oceanic_regrid.nc
#> names       : N2O.emissions 
#> zvar        : emi_n2o 

但是，正如所指出的，regrid输出（ out.raster）在CRS中显示为lon：我不确定为什么这样做，或者不确定 out.raster在范围上是全局的。

我还尝试通过两种方式限制剧情本身：

首先，我尝试在情节中添加 extents，即

plot(out.raster, ext=template.raster)
plot(map.us.proj, add=TRUE)

生成 the PDF的page = 2/4。不幸的是，这根本不会改变输出：AFAICS，它完全等于page = 1/4（上）。

其次，我尝试使用 raster::crop绑定重新定义的netCDF本身，然后使用与绑定重新定义的相同的 RasterLayer对象（ template.raster）进行绘制，然后再进行绘制：

template.raster.extent <- extent(template.raster)
out.raster.crop <-
  # overwrite the uncropped file
  crop(out.raster, template.raster.extent, filename=out.fp, overwrite=TRUE)
...
plot(out.raster.crop)
plot(map.us.proj, add=TRUE)

生成 the PDF的page = 3/4。不幸的是，它显然也完全等于page = 1/4（上）。

（如果您想知道，PDF的第4页是使用 fields::image.plot生成的，除非另有说明，否则该问题有一个不同的问题，描述为 here，除非StackOverflow重击了该链接。）

感谢您对这个R新手的帮助！

Answer 1

摘要

我的问题

How to limit display of an R `raster::plot` to the bounds of a `Raster*` object?

是基于我对问题的错误诊断。解决方案是正确设置底层 extent数据对象（即图的数据源）的边界（或 Raster*）


通过查询其基础文件来获取模板对象（用于创建数据对象）的边界（不做任何假设！）
通过查询其基础文件来获取模板对象的CRS（无需做任何假设！）
直接为模板对象设置边界和CRS


但


不要尝试仅在 Extent上设置分辨率；至少对我而言
raster::plot映射仍然存在一个较小的边界问题；至少相对于我在 fields::image.plot中使用的地图


细节

这个问题的答案实际上是另一个问题的答案：如何正确设置 extent对象的 Raster*？因为一旦我正确设置了要绘制的对象的范围， raster::plot问题就基本消失了。 （有一个小例外-见下文。）这可能是主要 raster开发人员Robert J. Hijmans试图用 this post传达的内容，但是如果是这样，我当时并没有意识到。

我得到了两个建议，这些建议本身并不是我所需要的，但它们使我走上了正确的道路，从而解决了这个问题。在这种情况下，该路径导致了 R package M3，这对于处理CMAQ和WRF输入和输出的数据很有用。这些建议是


使用 project.lonlat.to.M3进行重新网格化任务。
绘图问题可能与球形投影的生成模型（CMAQ）的假设有关。


第二个建议似乎是合理的，因为我知道我正在使用 +ellps=WGS84设置CRS。 （请参见 this repository中，特别是 regrid.global.to.AQMEII.r中带有注释的R。）因此，我决定对此进行研究。

我知道第一个建议只会遇到困难（因为它只能预测要点），但是为了确定起见，只看了 M3 doc。 ToC中的以下功能立即引起了我的注意：


get.proj.info.M3，它不仅从模板文件返回球形的PROJ.4字符串，而且还返回了我认为我需要提供的虚假的东西和北边的字符串：

# use package=M3 to get CRS from template file
out.crs <- get.proj.info.M3(template.in.fp)
cat(sprintf('out.crs=%s\n', out.crs)) # debugging
# out.crs=+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000

get.grid.info.M3，可以稍作努力，返回模板文件的范围/范围：

extents.info <- get.grid.info.M3(template.in.fp)
extents.xmin <- extents.info$x.orig
extents.xmax <- max(
  get.coord.for.dimension(
    file=template.in.fp, dimension="col", position="upper", units="m")$coords)
extents.ymin <- extents.info$y.orig
extents.ymax <- max(
  get.coord.for.dimension(
    file=template.in.fp, dimension="row", position="upper", units="m")$coords)
template.extents <-
  extent(extents.xmin, extents.xmax, extents.ymin, extents.ymax)

然后可以在模板 Raster*上设置这些界限

template.in.raster <- raster(template.in.fp, ...)
template.raster <- projectExtent(template.in.raster, crs=out.crs)
template.raster@extent <- template.extents

并使用模板重新定义输入 Raster*

out.raster <-
  projectRaster(
    # give a template with extents--fast, but gotta calculate extents
    from=in.raster, to=template.raster, crs=out.crs,
    # give a resolution instead of a template? no, that hangs
#    from=in.raster, res=grid.res, crs=out.crs,
    method='bilinear', overwrite=TRUE, format='CDF',
    # args from writeRaster
    NAflag=-999.0,  # match emi_n2o:missing_value,_FillValue (TODO: copy)
    varname=data.var.name, 
    varunit='ton N2O-N/yr',
    longname='N2O emissions',
    xname='COL',
    yname='ROW',
    filename=out.fp)
# above fails to set CRS, so
out.raster@crs <- CRS(out.crs)

（如上所述，通过设置模板进行重新栅格化仅用了7秒，但是在我杀死工作后的2小时后，通过设置网格分辨率进行的栅格化未能完成。）之后， raster::plot

map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
map.us.proj <-
  spTransform(map.us.unproj, CRS(out.crs)) # projected
...
pdf(file=pdf.fp, width=5.5, height=4.25)
...
plot(out.raster, # remaining args from image.plot
  main=title, sub=subtitle,
  xlab='', ylab='', axes=F, col=colors(100),
  axis.args=list(at=quantiles, labels=quantiles.formatted))
# add a projected CONUS map
plot(map.us.proj, add=TRUE)

几乎达到了预期的效果，但有一个例外：地图超出了数据范围，延伸到了北部和南部（尽管没有延伸到东部和西部），从而导致该图像与域的已发布图像不够相似，例如 this。有趣的是，当我用 fields::image.plot作图时

# see code in
# https://github.com/TomRoche/GEIA_to_netCDF/blob/master/plotLayersForTimestep.r
plot.raster(
  raster=out.raster,
  title=title,
  subtitle=subtitle,
  q.vec=probabilities.vec,
  colors,
  map.cmaq
)

我没有这个问题：。因此，至少在目前，我可能会使用 fields::image.plot进行显示。