как отобразить спроецированную карту на R::lattice::layerplot?

резюме: я могу отобразить карту долгота-широта на lattice::layerplot, и я могу отобразить карту конформной конической формы Ламберта (LCC) на spam::image. Как отобразить карту LCC на lattice::layerplot? Ниже приведены примеры того, как этого не следует делать — помощь в исправлении приветствуется (или даже просто в отладке).

детали:

Я успешно использовал решетчатую графику (через latticeExtra и rasterVis) для отображения непроецированных глобальных данных о долготе и широте, которые работают достаточно хорошо (хотя я определенно заинтригован ggplot/ggmap). В частности, я могу накладывать на эти графики карты, что очень важно для той работы, которой я занимаюсь. Однако в настоящее время я не могу использовать lattice::layerplot для некоторых региональных данных, спроецированных LCC: графики данных, но я не могу получить карту для наложения. Я могу сделать это более грубыми способами, но предпочел бы знать, как это сделать в lattice/rasterVis или подобном (например, ggplot/ggmap). Далее следуют два почти автономных примера... но если вы уже знаете, как это сделать, сообщите мне об этом, и я пропущу отладку. (Я сильно отстал от проекта.)

Данные netCDF, ozone_lcc.nc, поставляются с CRAN package=M3... за исключением того, что M3 поставляет его как

system.file("extdata/ozone_lcc.ncf", package="M3")

и это расширение файла (.ncf) в настоящее время вызывает проблемы для CRAN package=raster (см. эту публикацию). Вы можете либо переименовать этот файл (и поместить его в текущий рабочий каталог), либо загрузить только этот файл (270 КБ), или вы можете получить его из Архив M3 (но не забудьте переименовать его!).

Затем вы можете запустить любой из следующих примеров (при условии, что (IIRC) вы не используете Windows, на которой package=M3 не будет собираться (но ICBW)), изменив константы по желанию для вашей системы. Пример 1 создает карту типа, который, как я знаю (из предыдущего опыта), будет работать с raster в levelplot; однако в этом случае значения координат карты и данных/растра не совпадают. Пример 2 использует базовую графику в старом стиле и на самом деле отображает как данные, так и карту; к сожалению, я не знаю, как сделать карту, которую он генерирует, наложением на levelplot. Поскольку я хотел бы, чтобы этот код работал с кучей других кодов, использующих raster и 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, которое хочет sp::sp.lines, а тем более к тому, что хочет latticeExtra::layer. (Я достаточно новичок, чтобы найти документацию по решетке довольно непостижимой.) Кроме того, я бы предпочел не выполнять описанное выше преобразование IOAPI вручную (хотя я, конечно, предпочел бы сделать это, чем прыгать через обручи графики в старом стиле). ).