まとめ
R の表示をraster::plot
Raster* オブジェクトの境界に制限する方法は? 私が尋ねる理由:
私はRの初心者です
- 非投影 (lon-lat) グローバル空間データを ASCII からnetCDFに変換(解決済み)
- 「リグリッド」: つまり、データを領域に制限し、投影し、ダウンスケール (グリッドセルのサイズを縮小) します (ほとんど解決済み)。
残念なことに、科学的な作業では、主に視覚検査によるこのようなデータ変換の QA が行われます。上記のステップ 1 は良さそうです。しかし、ステップ 2 のほうが見栄えがよくなりましたが、本当はregriddedの境界 (または)のみをプロットしたいと考えています。私は、変換ステップを実行し、変換された出力を明らかに正しくプロットするパブリック git リポジトリ (ここ) でbash スクリプトによって駆動される R コードを持っています。ただし、リグリッディングの出力をプロットしようとする私の試みは、まったく正しくありません (出力の最初の 3 ページはこちらを参照してください)。extents
RasterLayer
raster::plot
直し方?
詳細
バックグラウンド
いくつかのデータ (地球規模の海洋排出インベントリ (EI)) を取得し、それを他のデータ (主に他の EI) と組み合わせて、モデルに入力する必要があります。このモデルは netCDF を消費したいと考えており、その netCDF を隣接する US (CONUS) よりもわずかに大きい空間ドメインで望んでいます。この画像の境界線がドメインの境界です。(領域の一部だけが海洋性であることに注意してください。) モデルはまた、netCDF が特定の方法 ( 12 km の解像度でのLCC ) を投影することを望んでいます。私はこれを2つの独立した問題として扱っています:
- グローバル EI をネイティブ ASCII 形式から netCDF に変換します
- netCDF をグローバル/非投影から縮小 (より細かい解像度) の投影サブドメインに「再グリッド化」します。
このパブリック git リポジトリにアーカイブしたコードを使用して、これらの問題を解決しようとしています。その git リポジトリのクローンを作成し、READMEGEIA_to_netCDF.sh
の「最初の例」で説明されているように bash ドライバーを構成/実行すると (特に packages= ,を使用して R が適切に構成されていると仮定すると)、netCDF が出力され、( を使用して)プロットされます。これはPDFに出力され、次のようになります。出力データの分布は正しいように見えるので、問題 1 は解決されたようです。raster
ncdf4
fields::image.plot
問題
問題 2 (別名READMEの「2 番目の例」) を解決するには、R package= が必要なようですraster
。エラーなしで実行され、その出力の PDF が表示されますregrid_GEIA_netCDF.sh
。は、 の最後の 4 つのコード ブロックに対応する 4 ページで構成されます。ここで関連するのは最初の 3 ページのみです (4 ページ目は使用しようとしており、より大きな問題があります)。regrid.global.to.AQMEII.r
GEIA_N2O_oceanic_regrid.pdf
regrid.global.to.AQMEII.r
fields::image.plot
Page=1/4 は、再グリッド化された単純なraster::plot
出力に加えて、次の CONUS マップの投影バージョンから得られwrld_simpl
ます。
plot(out.raster)
plot(map.us.proj, add=TRUE)
残念ながら、出力はグローバルに、またはほぼグローバルに見えます。しかし、出力が再グリッド化された目的のドメインは、はるかに小さくなっています。だから私は3つの質問があります(1つの主要な質問、2つのフォローアップ):
raster::plot
デフォルトで表示される画像は、引数で指定された範囲のみを表示しますか?RasterLayer
- もしそうなら、私のregridの何が問題なのですか? (詳細は以下)
- そうでない場合 (つまり、
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 自体をバインドするために使用しようとしました。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 初心者への支援に感謝します。