nsidc_polar_area.R

library(raadtools)

## Loading required package: raster

## Loading required package: sp

x <- readice("2017-09-20")

library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

p <- spex::polygonize(x) %>% 
  st_transform("+proj=laea +lat_0=-90 +datum=WGS84") %>% 
  mutate(area = st_area(geometry))

ice_area <- setValues(x, as.numeric(p$area)/1e6)
plot(ice_area, col = viridis::viridis(100))
contour(ice_area[[1]], add = TRUE)
contour(is.na(x), add = TRUE)
## where is line of true scale? 
##"+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"
contour(setValues(ice_area[[1]], 
                  rgdal::project(coordinates(ice_area), projection(ice_area), inv = TRUE)[,2]), 
        levels = -70, add = TRUE, col = "white", lwd = 2)