R：多边形数据的极地图投影

Question

我有什么：

• 北极和南极点
• 来自北极和南极各种地球物理实体的栅格数据

我想要什么：

立体或任何其他极地投影的地图，带有背景地图或海岸线，裁剪到点的范围。换句话说：像上面这样的地图，上面有我自己选择的底图。

到目前为止我做了什么：

我加载了所有数据（包括来自 naturalearthdata 的地表数据；参见 MWE）projected them into stereographic 并绘制了它。包含多边形数据的结果如下所示：

我的 MWE：

library(raster)
library(sf)
library(ggplot2)
library(rgdal)


# file load ---------------------------------------------------------------

# sea ice raster data
if (!file.exists("seaiceraster.tif")) {
  url = "https://seaice.uni-bremen.de/data/smos/tif/20100514_hvnorth_rfi_l1c.tif"
  download.file(url, destfile = 'seaiceraster.tif')
}
si.raster = raster::raster('seaiceraster.tif')


# land surface shapefile
if (!file.exists("110m-admin-0-countries")) {
  url_land = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_land.zip"
  download.file(url_land, destfile = "110m-admin-0-countries")
  unzip("110m-admin-0-countries")
}
world_shp = rgdal::readOGR("ne_10m_land.shp")

# points
p.data =  structure(
  list(
    Lat = c(
      73.0114126168676,70.325555278764,77.467797903163,
      58.6423827457304,66.3616310851294,59.2097857474643,
      75.3135274436283,60.1983078512275,72.6614399747201,
      61.1566678672946,73.0822309615673,55.7759666826898,
      75.1651656433833,69.0130753414173,62.3288262448589
    ),
    Lon = c(
      -59.9175490701543,-80.1900239630732,-40.4609968914928,
      -61.0914448815381,-60.0703668488408,-21.027205418284,
      -100.200463810276,-74.861777073788,-55.1093773178206,
      -29.4108649230234,-64.5878251008461,-36.5343322019187,
      -31.647365623387,-67.466355105829,-64.1162329769077
    )
  ),
  row.names = c(
    1911L, 592L,2110L,3552L,3426L,1524L,635L,4668L,
    3945L,2848L,3609L,36L,4262L,3967L,2725L
  ),
  class = "data.frame"
)

p = sf::st_as_sf(p.data, coords = c("Lon", "Lat"),
                 crs = "+init=epsg:4326")

# project -----------------------------------------------------------------

polar.crs = CRS("+init=epsg:3995")

si.raster.proj = projectRaster(si.raster, crs = polar.crs)
world_shp.proj = sp::spTransform(world_shp, polar.crs)
p.proj = sf::st_transform(p, polar.crs)

# preparation -------------------------------------------------------------

AG = ggplot2::fortify(world_shp.proj)

# make raster to data.frame
si.raster.df = si.raster.proj %>%
  raster::crop(., p.proj) %>%
  raster::rasterToPoints(., spatial = TRUE) %>%
  as.data.frame(.)

colnames(si.raster.df) = c("val", "x", "y")

# plot --------------------------------------------------------------------

ggplot() +
  # geom_polygon(data = AG, aes(long, lat, group = group)) + # un-comment to see
  geom_raster(data = si.raster.df, aes(x = x, y = y, fill = val)) +
  geom_sf(data = p.proj, color = "green", size = 3)

Answer 1

我已经稍微更改了您示例中的工作流程，为海冰数据添加了 stars 包，但我认为它应该可以满足您的需求。您需要调整裁剪大小以将其扩大一点，因为点 p 就在绘图区域的边缘。 st_buffer 可能会对此有所帮助。

我对所有对象都使用了 seaicebuffer.tif 文件中的 crs。

.tif 文件有一个我无法在我的计算机上轻松转换的 crs。它似乎可以使用米作为长度单位，并且可能是极地立体（变体 B）投影。不过，积分和世界数据似乎没有问题，这就是我一直使用它的原因。

library(raster)
library(sf)
library(ggplot2)
library(rgdal)
library(stars)

si <- stars::read_stars('seaiceraster.tif')

world_sf = rgdal::readOGR("ne_10m_land.shp") %>% 
  st_as_sf() %>%
  st_transform(st_crs(si))

# p <- ... same as example and then:
p <- st_transform(p, st_crs(si))

# get a bounding box for the points to crop si & world.
p_bbox <- st_bbox(p) %>% 
  st_as_sfc() %>% 
  st_as_sf() %>% 
  st_buffer(100000)

# crop si & world_sf to an area around the points (p)
world_cropped <- st_crop(world_sf, p_bbox)
si_cropped <- st_crop(si, p_bbox)

#Plot 
ggplot() + 
  geom_sf(data = world_cropped, 
          color = 'black', 
          fill = 'NA', 
          size = .2) + 
  geom_stars(data = si_cropped) + 
  geom_sf(data = p, color = 'red') + 
  scale_fill_continuous(na.value = 0)

南方 .tif 的丑陋黑客，星号读取为因素：

si <- stars::read_stars('20150324_hvsouth_rfi_l1c.tif', NA_value = 0 )
si$"20150324_hvsouth_rfi_l1c.tif" <- as.numeric(si$"20150324_hvsouth_rfi_l1c.tif")

ggplot() + geom_stars(data = si)