1

I would like to put a layer on my basemap that shows sea ice coverage levels for both the Beaufort and Chukchi Seas off Alaska for the last few years. Is this available anywhere? Would I have to make it myself? If there is no premade layer, where would the data be available?

  • 1
    The National Ice Center would be a good place to start. – Vince Jun 12 '15 at 19:40
  • When you say coverage I'm not sure if you mean thickness or just area. If there isn't such an exiting dataset and all you care about is area coverage, you may be looking at a remote sensing question - generating polygons based on satellite or other imagery. – Chris W Jun 12 '15 at 20:37
2

You can get daily data at 25km resolution from the NSIDC, here using R. I always prefer to go back to the rawest data possible, and derive my own metrics. It's just a straightforward binary format but not many systems read it with full understanding of the spatial and temporal metadata. The daily files are understanded by raster.

## daily file, 0.3 Mb 25km resolution
f <- "ftp://sidads.colorado.edu/pub/DATASETS/nsidc0081_nrt_nasateam_seaice/north/nt_20150603_f17_nrt_n.bin"
fbase <- basename(f)
if (!file.exists(fbase)) download.file(f, fbase, mode = "wb")

library(raster)

## we can get the data in with 
ice <- raster(fbase)

ice
class       : RasterLayer 
dimensions  : 448, 304, 136192  (nrow, ncol, ncell)
resolution  : 25000, 25000  (x, y)
extent      : -3837500, 3762500, -5362500, 5837500  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs 
data source : in memory
names       : layer 
values      : 0, 100  (min, max)
time        : 2015-06-03 

## optionally write out to GeoTIFF (but beware of how your GIS handles the Hughes 1980 ellipsoid)
## writeRaster(ice, filename = gsub("bin", "tif", fbase))

Make a plot with a map

plot(ice, axes = FALSE)
library(rgdal)
library(maptools); data(wrld_simpl)

wm <- spTransform(subset(wrld_simpl, coordinates(wrld_simpl)[,2] > 30),     CRS(projection(ice)))
plot(wm, add = TRUE)
llgridlines(wm, ndiscr = 50, side = "EN")
title(fbase)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.