4

まとめ

R の表示をraster::plotRaster* オブジェクトの境界に制限する方法は? 私が尋ねる理由:

私はRの初心者です

  1. 非投影 (lon-lat) グローバル空間データを ASCII からnetCDFに変換(解決済み)
  2. 「リグリッド」: つまり、データを領域に制限し、投影し、ダウンスケール (グリッドセルのサイズを縮小) します (ほとんど解決済み)。

残念なことに、科学的な作業では、主に視覚検査によるこのようなデータ変換の QA が行われます。上記のステップ 1 は良さそうです。しかし、ステップ 2 のほうが見栄えがよくなりましたが、本当はregriddedの境界 (または)のみをプロットしたいと考えています。私は、変換ステップを実行し、変換された出力を明らかに正しくプロットするパブリック git リポジトリ (ここ) でbash スクリプトによって駆動される R コードを持っています。ただし、リグリッディングの出力をプロットしようとする私の試みは、まったく正しくありません (出力の最初の 3 ページはこちらを参照してください)。extentsRasterLayerraster::plot

直し方?

詳細

バックグラウンド

いくつかのデータ (地球規模の海洋排出インベントリ (EI)) を取得し、それを他のデータ (主に他の EI) と組み合わせて、モデルに入力する必要があります。このモデルは netCDF を消費したいと考えており、その netCDF を隣接する US (CONUS) よりもわずかに大きい空間ドメインで望んでいます。この画像AQMEII-NA ドメインの境界線がドメインの境界です。(領域の一部だけが海洋性であることに注意してください。) モデルはまた、netCDF が特定の方法 ( 12 km の解像度でのLCC ) を投影することを望んでいます。私はこれを2つの独立した問題として扱っています:

  1. グローバル EI をネイティブ ASCII 形式から netCDF に変換します
  2. netCDF をグローバル/非投影から縮小 (より細かい解像度) の投影サブドメインに「再グリッド化」します。

このパブリック git リポジトリにアーカイブしたコードを使用して、これらの問題を解決しようとしています。その git リポジトリのクローンを作成し、READMEGEIA_to_netCDF.shの「最初の例」で説明されているように bash ドライバーを構成/実行すると (特に packages= ,を使用して R が適切に構成されていると仮定すると)、netCDF が出力され、( を使用して)プロットされます。これはPDFに出力され、次のようになります。出力データの分布は正しいように見えるので、問題 1 は解決されたようです。rasterncdf4fields::image.plotGEIA 世界の海洋排出量

問題

問題 2 (別名READMEの「2 番目の例」) を解決するには、R package= が必要なようですraster。エラーなしで実行され、その出力の PDF が表示されますregrid_GEIA_netCDF.sh。は、 の最後の 4 つのコード ブロックに対応する 4 ページで構成されます。ここで関連するのは最初の 3 ページのみです (4 ページ目は使用しようとしており、より大きな問題があります)。regrid.global.to.AQMEII.rGEIA_N2O_oceanic_regrid.pdfregrid.global.to.AQMEII.rfields::image.plot

Page=1/4 は、再グリッド化された単純なraster::plot出力に加えて、次の CONUS マップの投影バージョンから得られwrld_simplます。

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

残念ながら、出力はグローバルに、またはほぼグローバルに見えます。raster::プロット境界の問題しかし、出力が再グリッド化された目的のドメインは、はるかに小さくなっていますAQMEII-NA ドメイン。だから私は3つの質問があります(1つの主要な質問、2つのフォローアップ):

  1. raster::plotデフォルトで表示される画像は、引数で指定された範囲のみを表示しますか?RasterLayer
  2. もしそうなら、私のregridの何が問題なのですか? (詳細は以下)
  3. そうでない場合 (つまり、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 は、ここに示す出力ドメインの定義と一致しているように見えることに注意してください。(リンク先のページに移動したら、地図をスクロールします。)

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 の出力 ( ) は、CRS で lon-lat であることを示しています。なぜそうなのか、またはエクステントがグローバルでout.rasterあることを意味するのかはわかりません。out.raster

また、次の 2 つの方法でプロット自体を制限しようとしました。

まず、プロットに を追加してみましextentsた。つまり、

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

PDFの page=2/4 を生成します。残念ながら、それは出力をまったく変更しません。AFAICS、これは page=1/4 (上記) と完全に同じです。

次に、リグリッドのバインドに使用したのと同じオブジェクト ( ) を使用raster::cropして、それをプロットする前に、再グリッド化された netCDF 自体をバインドするために使用しようとしました。RasterLayertemplate.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;に解像度を設定しようとしないでください。少なくとも私にとっては、それはハングアップします
  • マップにはまだ1つの小さな境界の問題がありraster::plotます。少なくとも、私が使用しているマップと比較してfields::image.plot

詳細

この質問への答えは、実際には別の質問への答えです。オブジェクトのを正しく設定extentする方法は?プロットしたいオブジェクトの範囲を適切に設定するとRaster*、問題はほとんど解消されたからです。raster::plot(1つの小さな例外を除いて-以下を参照してください。)これは、メインの開発者であるRobert J. Hijmansがこの投稿rasterで伝えようとしていたことかもしれませんが、そうだとすれば、当時私はそれを認識していませんでした。

私は、それ自体は必要なものではなかったが、適切な道に私を導いた2つの提案を得た後、問題を解決しました。この場合、そのパスはRパッケージM3につながりました。これは、CMAQおよびWRFとの間で送受信されるデータを処理するのに役立ちます。提案は

  1. project.lonlat.to.M3再グリッドタスクに使用します。
  2. プロットの問題は、球面投影の生成モデル(CMAQ)による仮定に関連している可能性があります。

2番目の提案は、私がCRSを取得して設定していることを知っていたので、もっともらしいように思われました+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)

1つの例外を除いて、ほぼ予想どおりでした。ラスター::プロットを使用した過剰な南北プロット範囲マップはデータの境界を超えて北と南に広がり(東と西ではありません)、画像がドメインの公開された画像に十分に似ていない原因になっています。興味深いことに、私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。ですからfields::image.plot、少なくとも当面は、おそらくディスプレイに使用します。

于 2012-11-01T00:40:44.267 に答える