R: Polar map projection of polygon data

Question

What I have:

• points in the arctic and antarctic
• raster data from various geophysical entities in arctic and antarctic

What I want:

A map in stereographic or any other polar projection with background map or coastlines, cropped to the extent of the points. In other words: A map like above with base map of my own choice.

What I did so far:

I loaded all the data (including land surface data from naturalearthdata; see MWE), projected them into stereographic and plotted that. The result including the polygon data looks then like this:

My 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

I've changed the workflow in your example a bit to add the stars package for the sea ice data, but I think it should get you what you're looking for. You'll need to adjust the crop size to expand it a little, as the points p are right on the edge of the plotted area. st_buffer might help with that.

I used the crs from the seaicebuffer.tif file for all of the objects.

The .tif file has a crs that I'm not able to easily transform on my computer. It seems to be able to use meters as a lengthunit and might be a polar stereographic (variant B) projection. The points & world data don't seem to have a problem transforming to it though, which is why I've used it throughout.

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)

Ugly hack for the southern .tif that stars reads as factors:

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)