knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "README-" )
The goal of sospatial is to provide some easy to use Southern Ocean data sets for mapping and exploration in R.
We aim to fill some gaps that make mapping in polar regions more difficult that usual.
This example shows the used of the packed_lats
cache to draw a rough map of the intensity
of where the daily sea ice extent was in the Southern Ocean.
WIP - when this data set is fully spatially and temporally referenced we can use it to extract abstract metrics like "max sea ice extent", which might mean the major 'edge' on the day of greatest extent, or the cumulatively most-northerly location at all longitudes during a year ...
## reconstruct record library(dplyr) library(sospatial) ice <- tibble(lon = rep_len(seq(-180, 179), length.out = length(packed_lats)), lat = packed_lats/10, day = rep(ice_dates, each = 360)) prjstere <- "+proj=stere +lat_0=-90 +datum=WGS84" ice[c("x", "y")] <- rgdal::project(as.matrix(ice[c("lon", "lat")]), prjstere) library(ggplot2) ggplot(ice, aes(x, y)) + geom_bin2d(bins = 120) + coord_equal() # what does a month look like mon <- ice %>% dplyr::filter(format(day, "%m") == "08") ggplot(mon, aes(x, y)) + geom_bin2d(bins = 120) + coord_equal() + xlim(range(mon$x)) + ylim(range(mon$y)) mon$era <- c("old", "new")[(mon$day > as.POSIXct("1998-06-15")) + 1] ggplot(mon, aes(x, y)) + geom_bin2d(bins = 120) + coord_equal() + xlim(range(mon$x)) + ylim(range(mon$y)) + facet_wrap(~era) ## add the fronts library(orsifronts) library(dplyr) fronts <- ggplot2::fortify(sp::spTransform(orsifronts, prjstere)) ggplot(mon, aes(x, y)) + geom_bin2d(bins = 120) + geom_path(data = fronts, aes(long, lat, group = group, colour = id)) ## get a coastline cst <- fortify(sp::spTransform(rnaturalearth::ne_coastline(), prjstere)) ggplot(mon, aes(x, y)) + geom_bin2d(bins = 120) + geom_path(data = fronts, aes(long, lat, group = group, colour = id)) + geom_path(data = cst, aes(long, lat, group = group)) + xlim(range(fronts$long)) + ylim(range(fronts$lat))
Pull out a specific summary line, say the maximum latitude per longitude for any day in October in the last ten years.
This is a bit of a dummy example, but is the kind of task that we hope to make easier so that these basic summaries can be produced by anyone.
oct <- ice %>% dplyr::filter(day > as.POSIXct("2007-01-01"), format(day, "%m") == "10") %>% filter(!is.na(lat)) %>% group_by(lon) %>% summarize(lat = max(lat)) ## construct that as a line, and a polygon library(sf) line <- st_sf(geometry = st_sfc(st_linestring(cbind(oct$lon, oct$lat)), crs = 4326), name = "october_max_2007_2017") poly <- st_sf(geometry = st_polygonize(st_transform(st_sfc(st_linestring(cbind(oct$lon, oct$lat)[c(1:nrow(oct), 1), ]), crs = 4326), prjstere)), name = "october_max_2007_2017") plot(st_geometry(line)) plot(orsifronts, add = TRUE, col = "firebrick") plot(st_geometry(poly)) ## area is a placeholder here, obviously we need to cut out the continent ## and any polynas ## area in stere st_area(poly) ## (more) true area ## note that this is more like what st_area will give for a longlat dataset ## since stereographic is not area-preserving and laea will be more faithful ## the ellipsoidal methods used for longlat st_area(st_transform(poly, "+proj=laea +lat_0=-90 +datum=WGS84"))
Please note that this project is released with a Contributor Code of Conduct. By participating in this project you agree to abide by its terms.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.