如何在R :: lattice :: layerplot上显示投影地图? [英] how to display a projected map on an R::lattice::layerplot?
问题描述
摘要:我可以在lattice::layerplot上显示lon-lat映射，并且可以在spam::image上显示Lambert保形圆锥(LCC)映射.如何在lattice::layerplot上显示LCC映射?以下是如何不执行此操作的示例-赞赏了修复方面的帮助(甚至只是调试).
summary: I can display a lon-lat map on a lattice::layerplot, and I can display a Lambert conformal conic (LCC) map on a spam::image. How to display an LCC map on a lattice::layerplot? Examples of how not to do it follow--assistance in fixing is appreciated (or even just debugging).
详细信息:
我一直在成功使用网格图形(通过latticeExtra和rasterVis)来显示未投影的经纬度全球大气数据，效果很好(尽管我对ggplot/ggmap确实很感兴趣) .特别是，我能够用地图覆盖这些图，这对于我正在做的工作是必不可少的.但是，我目前无法对某些区域数据投影的LCC使用lattice::layerplot:数据图，但是我无法覆盖地图.我可以通过更粗略的方式来执行此操作，但希望了解如何使用lattice/rasterVis或类似的方式(例如ggplot/ggmap)执行此操作.接下来是两个几乎独立的示例...但是，如果您已经知道如何执行此操作，请告诉我，我将跳过调试. (我在一个项目上落后了.)
I've been using trellis graphics (via latticeExtra and rasterVis) successfully to display unprojected lon-lat global atmospheric data, which works well enough (though I am definitely intrigued by ggplot/ggmap). Particularly, I am able to overlay those plots with maps, which is essential for the sort of work I'm doing. However I'm currently unable to use lattice::layerplot for some regional data projected LCC: the data plots, but I cannot get a map to overlay. I can do this by cruder means, but would prefer to know how to do this in lattice/rasterVis or similar (e.g., ggplot/ggmap). Two nearly-self-contained examples follow ... but if you know how to do this already, please let me know, and I'll skip debugging. (I'm waayyy behind on a project.)
netCDF数据ozone_lcc.nc随 CRAN package = M3 ...，除了M3将其提供为
The netCDF data, ozone_lcc.nc, comes with CRAN package=M3 ... except that M3 supplies it as
system.file("extdata/ozone_lcc.ncf", package="M3")
，并且该文件扩展名(.ncf)当前导致 CRAN package =光栅(请参见这篇文章).您可以重命名该文件(并将其放在当前工作目录中)，也可以下载只是该文件(270 kB)，也可以从
and that file extension (.ncf) currently causes problems for CRAN package=raster (see this post). You can either rename that file (and put it some current working directory), or you can download just that file (270 kB), or you can get it from the M3 tarball (but remember to rename it!).
然后，您可以运行以下任何示例(未运行Windows的情况下提供的(IIRC)，在Windows上不会生成package = M3(但为ICBW))，可以根据需要更改常量以适应您的系统.示例1生成了一种类型的映射，我知道(根据以前的经验)，该映射将与levelplot中的raster一起使用；但是，在这种情况下，地图的坐标值与数据/栅格的坐标值不匹配.示例2使用了较旧样式的基础图形，并同时绘制了数据和地图.不幸的是，我不知道如何制作在levelplot上生成的地图.正如我想使此代码与其他使用raster和levelplot的代码一起使用时，这是一个问题.
You can then run any of the following examples (provided (IIRC) you're not running Windows, on which package=M3 won't build (but ICBW)), altering constants as desired to accommodate your system. Example 1 generates a map of a type that I know (from previous experience) will work with a raster in a levelplot; however, in this case, the coordinate values of the map and the data/raster don't match. Example 2 uses older-style base graphics, and actually plots both data and map; unfortunately, I don't know how to make the sort of map which it generates overlay on a levelplot. As I'd like to make this code work with a bunch of other codes which use raster and levelplot, that's a problem.
示例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返回的简单矩阵(等，但不多)到sp::sp.lines想要的SpatialLines，更不用说latticeExtra::layer想要的东西了. (我是一个新手，足以找到难以理解的格栅文档.)此外，我宁愿避免手动进行上述IOAPI转换(尽管我当然愿意这样做，而不是跳过旧式图形的框架) ).
But I can't see how to get from the simple matrix (et al, but not much) returned by M3::get.map.lines.M3.proj to the SpatialLines that sp::sp.lines wants, much less to whatever latticeExtra::layer wants. (I'm enough of a newbie to find the trellis docs rather inscrutable.) Furthermore I'd prefer to avoid doing the IOAPI conversion above manually (though I'd certainly rather do that than jump through the hoops of the old-style graphics).
推荐答案
做到这一点的一种方法虽然很丑，但是却是使用线性IOAPI变换将"maps::map"中的坐标固定"在"c22"中.例如，
One way to do this, though ugly, is to "fix" the coordinates in the state.map obtained from maps::map, using the linear IOAPI transform. E.g.,
示例3:产生类似
##### start example 3 #####
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")
# See notes in question regarding unfortunate problem with raster::raster,
# and remember to download or rename the file (symlinking alone will not work).
lcc.file <- "./ozone_lcc.nc"
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.example3.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
# replace its coordinates with
metadata.coords.IOAPI.list <- M3::get.grid.info.M3(lcc.file)
metadata.coords.IOAPI.x.orig <- metadata.coords.IOAPI.list$x.orig
metadata.coords.IOAPI.y.orig <- metadata.coords.IOAPI.list$y.orig
metadata.coords.IOAPI.x.cell.width <- metadata.coords.IOAPI.list$x.cell.width
metadata.coords.IOAPI.y.cell.width <- metadata.coords.IOAPI.list$y.cell.width
map.lines <- M3::get.map.lines.M3.proj(
file=lcc.file, database="state", units="m")
map.lines.coords.IOAPI.x <-
(map.lines$coords[,1] - metadata.coords.IOAPI.x.orig)/metadata.coords.IOAPI.x.cell.width
map.lines.coords.IOAPI.y <-
(map.lines$coords[,2] - metadata.coords.IOAPI.y.orig)/metadata.coords.IOAPI.y.cell.width
map.lines.coords.IOAPI <-
cbind(map.lines.coords.IOAPI.x, map.lines.coords.IOAPI.y)
# start debugging
class(map.lines.coords.IOAPI)
summary(map.lines.coords.IOAPI)
# end debugging
# > class(map.lines.coords.IOAPI)
# [1] "matrix"
# > summary(map.lines.coords.IOAPI)
# map.lines.coords.IOAPI.x map.lines.coords.IOAPI.y
# Min. : 12.88 Min. :14.47
# 1st Qu.: 78.62 1st Qu.:39.28
# Median :101.37 Median :57.25
# Mean : 95.17 Mean :55.65
# 3rd Qu.:124.47 3rd Qu.:72.57
# Max. :140.51 Max. :93.16
# NA's :168 NA's :168
coordinates(state.map.shp) <- map.lines.coords.IOAPI # FAIL:
> Error in `coordinates<-`(`*tmp*`, value = c(99.0242231482288, 99.1277727928691, :
> setting coordinates cannot be done on Spatial objects, where they have already been set
state.map.IOAPI <- state.map # copy
state.map.IOAPI$x <- map.lines.coords.IOAPI.x
state.map.IOAPI$y <- map.lines.coords.IOAPI.y
state.map.IOAPI$range <- c(
min(map.lines.coords.IOAPI.x),
max(map.lines.coords.IOAPI.x),
min(map.lines.coords.IOAPI.y),
max(map.lines.coords.IOAPI.y))
state.map.IOAPI.shp <-
maptools::map2SpatialLines(state.map.IOAPI, proj4string=lcc.crs)
# start debugging
# thanks, Felix Andrews!
class(state.map.IOAPI.shp)
summary(do.call("rbind",
unlist(coordinates(state.map.IOAPI.shp), recursive=FALSE)))
# end debugging
# > class(state.map.IOAPI.shp)
# [1] "SpatialLines"
# attr(,"package")
# [1] "sp"
# > summary(do.call("rbind",
# + unlist(coordinates(state.map.IOAPI.shp), recursive=FALSE)))
# V1 V2
# Min. : 12.88 Min. :14.47
# 1st Qu.: 78.62 1st Qu.:39.28
# Median :101.37 Median :57.25
# Mean : 95.17 Mean :55.65
# 3rd Qu.:124.47 3rd Qu.:72.57
# Max. :140.51 Max. :93.16
pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
sp::sp.lines(state.map.IOAPI.shp, lwd=0.8, col='darkgray'))
dev.off()
# change this for viewing PDF on your system
system(sprintf("xpdf %s", lcc.pdf))
##### end example 3 #####
这篇关于如何在R :: lattice :: layerplot上显示投影地图?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
