概要:に経度緯度マップをlattice::layerplot表示でき、 にランベルト正角円錐 (LCC) マップを表示できますspam::image。にLCCマップを表示するにはlattice::layerplot？やってはいけないことの例を以下に示します - 修正の支援を歓迎します (または単にデバッグするだけでも)。
詳細：
私はトレリス グラフィックスを ( および を介して) 使用しlatticeExtraてrasterVis、投影されていない lon-lat グローバル大気データを正常に表示しました。特に、これらのプロットをマップに重ねることができます。これは、私が行っている種類の作業に不可欠です。ただし、現在、LCC を投影した一部の地域データには使用できません。データ プロットですが、オーバーレイするマップを取得できません。これはもっと大雑把な方法で行うこともできますが、または同様の方法 (例: )で行う方法を知りたいと考えています。ほぼ自己完結型の例を 2 つ以下に示します...しかし、これを行う方法を既に知っている場合は、お知らせください。デバッグは省略します。（私はプロジェクトにかなり遅れています。）ggplotggmaplattice::layerplotlattice/rasterVisggplot/ggmap
netCDF データ にはCRAN package=M3ozone_lcc.ncが付属しています ... ただし、M3 はそれを次のように提供します。
system.file("extdata/ozone_lcc.ncf", package="M3")
そのファイル拡張子 ( ) は現在、 CRAN package=raster.ncfで問題を引き起こしています (この投稿を参照してください)。そのファイルの名前を変更する (そして現在の作業ディレクトリに置く) か、そのファイル(270 kB) だけをダウンロードするか、M3 tarballから取得することができます(ただし、名前を変更することを忘れないでください!)。
次に、次の例のいずれかを実行できます (ただし、(IIRC) Windows を実行していない場合、package=M3はビルドされません (ただし ICBW))、必要に応じてシステムに合わせて定数を変更します。例 1 は、(以前の経験から) araster内の a で動作することがわかっている型のマップを生成しlevelplotます。ただし、この場合、マップとデータ/ラスターの座標値は一致しません。例 2 は古いスタイルのベース グラフィックを使用し、実際にはデータとマップの両方をプロットします。残念ながら、私はそれが生成するような地図をlevelplot. rasterこのコードを、 andを使用する他の多くのコードで動作させたいのでlevelplot、それは問題です。
例 1: 次のような出力を生成します
##### start example 1 #####
library("M3") # http://cran.r-project.org/web/packages/M3/
library("rasterVis") # http://cran.r-project.org/web/packages/rasterVis/
## Use an example file with projection=Lambert conformal conic.
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3")
lcc.file <- "./ozone_lcc.nc" # unfortunate problem with raster::raster
lcc.proj4 <- M3::get.proj.info.M3(lcc.file)
lcc.proj4 # debugging
# [1] "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000"
# Note +lat_0=40 +lat_1=33 +lat_2=45 for maps::map@projection (below)
lcc.crs <- sp::CRS(lcc.proj4)
lcc.pdf <- "./ozone_lcc.example1.pdf" # for output
## Read in variable with daily ozone.
# o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs)
# ozone_lcc.nc has 4 timesteps, choose 1 at random
o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs, level=1)
o3.raster@crs <- lcc.crs # why does the above assignment not take?
# start debugging
o3.raster
summary(coordinates(o3.raster)) # thanks, Felix Andrews!
M3::get.grid.info.M3(lcc.file)
# end debugging
# > o3.raster
# class : RasterLayer
# band : 1
# dimensions : 112, 148, 16576 (nrow, ncol, ncell)
# resolution : 1, 1 (x, y)
# extent : 0.5, 148.5, 0.5, 112.5 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000
# data source : .../ozone_lcc.nc
# names : O3
# z-value : 1
# zvar : O3
# level : 1
# > summary(coordinates(o3.raster))
# x y
# Min. : 1.00 Min. : 1.00
# 1st Qu.: 37.75 1st Qu.: 28.75
# Median : 74.50 Median : 56.50
# Mean : 74.50 Mean : 56.50
# 3rd Qu.:111.25 3rd Qu.: 84.25
# Max. :148.00 Max. :112.00
# > M3::get.grid.info.M3(lcc.file)
# $x.orig
# [1] -2736000
# $y.orig
# [1] -2088000
# $x.cell.width
# [1] 36000
# $y.cell.width
# [1] 36000
# $hz.units
# [1] "m"
# $ncols
# [1] 148
# $nrows
# [1] 112
# $nlays
# [1] 1
# The grid (`coordinates(o3.raster)`) is integers, because this
# data is stored using the IOAPI convention. IOAPI makes the grid
# Cartesian by using an (approximately) equiareal projection (LCC)
# and abstracting grid positioning to metadata stored in netCDF
# global attributes. Some of this spatial metadata can be accessed
# by `M3::get.grid.info.M3` (above), and some can be accessed via
# the coordinate reference system (`CRS`, see `lcc.proj4` above)
## Get US state boundaries in projection units.
state.map <- maps::map(
database="state", projection="lambert", par=c(33,45), plot=FALSE)
# parameters to lambert: ^^^^^^^^^^^^
# see mapproj::mapproject
state.map.shp <-
maptools::map2SpatialLines(state.map, proj4string=lcc.crs)
# start debugging
# thanks, Felix Andrews!
class(state.map.shp)
summary(do.call("rbind",
unlist(coordinates(state.map.shp), recursive=FALSE)))
# end debugging
# > class(state.map.shp)
# [1] "SpatialLines"
# attr(,"package")
# [1] "sp"
# > summary(do.call("rbind",
# + unlist(coordinates(state.map.shp), recursive=FALSE)))
# V1 V2
# Min. :-0.234093 Min. :-0.9169
# 1st Qu.:-0.000333 1st Qu.:-0.8289
# Median : 0.080378 Median :-0.7660
# Mean : 0.058492 Mean :-0.7711
# 3rd Qu.: 0.162993 3rd Qu.:-0.7116
# Max. : 0.221294 Max. :-0.6343
# I can't see how to relate these coordinates to `coordinates(o3.raster)`
pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
dev.off()
# change this for viewing PDF on your system
system(sprintf("xpdf %s", lcc.pdf))
# data looks correct, map invisible
## Try again, with lambert(40,33)
state.map <- maps::map(
database="state", projection="lambert", par=c(40,33), plot=FALSE)
# parameters to lambert: ^^^^^^^^^^^^
# see mapproj::mapproject
state.map.shp <-
maptools::map2SpatialLines(state.map, proj4string=lcc.crs)
# start debugging
summary(do.call("rbind",
unlist(coordinates(state.map.shp), recursive=FALSE)))
# end debugging
# > summary(do.call("rbind",
# + unlist(coordinates(state.map.shp), recursive=FALSE)))
# V1 V2
# Min. :-0.2226344 Min. :-0.9124
# 1st Qu.:-0.0003149 1st Qu.:-0.8295
# Median : 0.0763322 Median :-0.7706
# Mean : 0.0553948 Mean :-0.7752
# 3rd Qu.: 0.1546465 3rd Qu.:-0.7190
# Max. : 0.2112617 Max. :-0.6458
# no real change from previous `coordinates(state.map.shp)`
pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
dev.off()
system(sprintf("xpdf %s", lcc.pdf))
# as expected, same as before: data looks correct, map invisible
##### end example 1 #####
例 2: 次のような出力を生成します
##### start example 2 #####
# Following adapted from what is installed in my
# .../R/x86_64-pc-linux-gnu-library/2.14/m3AqfigExampleScript.r
# (probably by my sysadmin), which also greatly resembles
# https://wiki.epa.gov/amad/index.php/R_packages_developed_in_AMAD
# which is behind a firewall :-(
## EXAMPLE WITH LAMBERT CONIC CONFORMAL FILE.
library("M3")
library("aqfig") # http://cran.r-project.org/web/packages/aqfig/
## Use an example file with LCC projection: either local download ...
lcc.file <- "./ozone_lcc.nc"
## ... or as installed by package=M3:
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3")
## Choose the one that works for you.
lcc.pdf <- "./ozone_lcc.example2.pdf" # for output
## READ AND PLOT OZONE FROM FILE WITH LCC PROJECTION.
## Read in variable with daily ozone. Note that we can give dates
## rather than date-times, and that will include all time steps
## anytime during those days. Or, we can give lower and upper bounds
## and all time steps between these will be taken.
o3 <- M3::get.M3.var(
file=lcc.file, var="O3", hz.units="m",
ldatetime=as.Date("2001-07-01"), udatetime=as.Date("2001-07-04"))
# start debugging
class(o3)
summary(o3)
summary(o3$x.cell.ctr)
# end debugging
# > class(o3)
# [1] "list"
# > summary(o3)
# Length Class Mode
# data 66304 -none- numeric
# data.units 1 -none- character
# x.cell.ctr 148 -none- numeric
# y.cell.ctr 112 -none- numeric
# hz.units 1 -none- character
# rows 112 -none- numeric
# cols 148 -none- numeric
# layers 1 -none- numeric
# datetime 4 POSIXct numeric
# > summary(o3$x.cell.ctr)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# -2718000 -1395000 -72000 -72000 1251000 2574000
# Note how these grid coordinates relate to the IOAPI metadata above:
# min(o3$x.cell.ctr) == -2718000
# == -2736000 + (36000/2) == x.orig + (x.cell.width/2)
## Get colors and map boundaries for plot.
library("fields")
col.rng <- tim.colors(20)
detach("package:fields")
## Get state boundaries in projection units.
map.lines <- M3::get.map.lines.M3.proj(
file=lcc.file, database="state", units="m")
# start debugging
class(map.lines)
summary(map.lines)
summary(map.lines$coords)
# end debugging
# > class(map.lines)
# [1] "list"
# > summary(map.lines)
# Length Class Mode
# coords 23374 -none- numeric
# units 1 -none- character
# > summary(map.lines$coords)
# x y
# Min. :-2272238 Min. :-1567156
# 1st Qu.: 94244 1st Qu.: -673953
# Median : 913204 Median : -26948
# Mean : 689997 Mean : -84644
# 3rd Qu.: 1744969 3rd Qu.: 524531
# Max. : 2322260 Max. : 1265778
# NA's :168 NA's :168
## Set color boundaries to encompass the complete data range.
z.rng <- range(as.vector(o3$data))
## Make a color tile plot of the ozone for the first day (2001-07-01).
pdf(file=lcc.pdf)
image(o3$x.cell.ctr, o3$y.cell.ctr, o3$data[,,1,1],
col=col.rng, zlim=z.rng,
xlab="x-coord (m)", ylab="y-coord (m)")
## Put date-time string and chemical name (O3) into a format I can use
## to label the actual figure.
date.str <- format(o3$datetime[1], "%Y-%m-%d")
title(main=bquote(paste(O[3], " on ", .(date.str), sep="")))
## Put the state boundaries on the plot.
lines(map.lines$coords)
## Add legend to right of plot. NOTE: YOU CANNOT ADD TO THE PLOT
## AFTER USING vertical.image.legend(). Before making a new plot,
## open a new device or turn this device off.
vertical.image.legend(zlim=z.rng, col=col.rng)
dev.off() # close the plot if you haven't already, ...
# ... and change the following to display PDFs on your system.
system(sprintf("xpdf %s", lcc.pdf))
# data displays with state map
##### end example 2 #####
しかし、単純な行列 (その他、しかしそれほど多くはありません) から返されM3::get.map.lines.M3.proj たものSpatialLinesから取得する方法がわかりません。(私はトレリスのドキュメントがかなり不可解であることに気付くのに十分な初心者です。) さらに、上記の IOAPI 変換を手動で行うことは避けたいと思います (ただし、古いスタイルのグラフィックスのフープを飛び越えるよりも、そうする方が確実です)。 ）。sp::sp.lineslatticeExtra::layer
