4

概括

如何将 Rraster::plot的显示限制在 Raster* 对象的范围内？为什么我问：

我是 R 初学者，必须

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

不幸的是，在科学工作中，主要是通过视觉检查来进行诸如此类的 QAs 数据转换。上面的第 1 步看起来不错。但是尽管第 2 步现在看起来好多了，但我真的只想绘制regridded的边界（或extents）RasterLayer。我在公共 git 存储库（此处）中有由 bash 脚本驱动的 R 代码，它执行转换步骤，并绘制转换后的输出，显然是正确的。但是，我尝试绘制 regridding 的输出raster::plot并不完全正确（请参阅此处的前 3 页输出）。

怎么修？

细节

背景

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

  1. 将全局 EI 从其本机 ASCII 格式转换为 netCDF
  2. 将 netCDF 从全局/未投影到缩小（更精细）的投影子域“重新网格化”。

我正在尝试使用我在这个公共 git 存储库中存档的代码来解决这些问题。如果您克隆该 git 存储库，然后按照 READMEGEIA_to_netCDF.sh中“第一个示例”的描述配置/运行 bash 驱动程序（并假设您的 R 已正确配置，特别是使用 packages= , ），它将输出 netCDF 并绘图（使用）输出到PDF，应该如下所示：输出数据的分布似乎是正确的，因此问题 1 似乎已解决。rasterncdf4fields::image.plotGEIA 全球海洋排放量

问题

解决问题 2（又名README 中的“第二个示例” ）似乎需要 R package= raster。我的 bash 驱动程序regrid_GEIA_netCDF.sh调用regrid.global.to.AQMEII.r（都在存储库中），它运行没有错误，并显示其输出的 PDF。GEIA_N2O_oceanic_regrid.pdf包含 4 页，对应 . 末尾的 4 个代码块regrid.global.to.AQMEII.r。这里只有前3页相关（第4页尝试使用fields::image.plot并且有更大的问题）。

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

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

不幸的是，输出看起来是全局的，或者几乎是全局的：raster::plot 边界问题但是输出重新网格化到的所需域要小得多：AQMEII-NA结构域. 所以我有 3 个问题（1 个主要问题，2 个后续问题）：

  1. raster::plot默认情况下显示的图像是否仅显示其参数中给定的范围RasterLayer？
  2. 如果是这样，我的 regrid 有什么问题？（更多下文）
  3. 如果不是（即，如果raster::plot默认显示超过其范围RasterLayer），如何使raster::plot仅显示范围RasterLayer？

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

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

请注意，CRS 似乎与此处给出的输出域的定义相匹配。（一旦在链接页面，滚动过去的地图。）

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

但是，如前所述，regridout.raster输出out.raster（

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

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

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

生成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)

生成PDF的 page=3/4 。不幸的是，它显然也与 page=1/4（上图）完全相同。

（如果您想知道，PDF 的第 4 页是使用 生成的fields::image.plot，它有一个不同的问题，在此处描述，除非 StackOverflow 破坏了该链接。）

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

4

1 回答 1

3

概括

我的问题

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

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

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

但

  • 不要试图只在Extent;上设置分辨率 至少对我来说，那挂起
  • 地图还有一个小的边界问题raster::plot；至少，相对于我正在使用的地图fields::image.plot

细节

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

在得到两个建议后，我解决了这个问题，这些建议本身并不是我需要的，但它们让我走上了正确的道路。在这种情况下，该路径通向R 包M3，这对于处理 CMAQ 和 WRF 的数据输入和输出非常有用。这些建议是

  1. 用于project.lonlat.to.M3重新网格化任务。
  2. 绘图问题可能与球面投影的生成模型 (CMAQ) 的假设有关。

第二个建议似乎是合理的，因为我知道我正在使用+ellps=WGS84. （请参阅此存储库中大量评论的 R ，特别是regrid.global.to.AQMEII.r。）所以我决定研究一下。

我知道第一个建议只会有困难（因为它只投射点数），但为了确定，我查看了M3 文档。它的 ToC 中的以下功能立即引起了我的注意：

  1. 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

  2. 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)

几乎与预期的一样，但有一个例外：raster::plot 的南北绘图范围过多地图向北和向南（尽管不是向东和向西）延伸到数据的边界之外，导致图像与该域的已发布图像不够相似，例如this。有趣的是，当我用fields::image.plotlike

# 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。所以我可能会用它fields::image.plot来展示，至少目前是这样。

于 2012-11-01T00:40:44.267 回答