I want to plot remote sensed sea ice concentration data for Antarctica, in Program R. The product website (http://cersat.ifremer.fr/data/products/catalogue) says the data are in Antarctic polar stereographic projection. However, the extent of the rasters (thus the x/y (lat/lon) coordinates) does not agree with this projection.
Preamble. This works, and shows what I expect:
This is a bathymetry map of the Antarctic Peninsula, in WGS84 projection.
I can reproject this map to Antarctic polar stereographic projection. In R:
Polar.prj = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs "
Site.polar = projectRaster(from = site.WGS84, crs = Polar.prj)
I end up with a raster with extent values in the “negative millions” (longitude) and “positive millions” (latitude) (see map image, and raster details below)
> Site.polar
class : RasterLayer
dimensions : 1321, 3217, 4249657 (nrow, ncol, ncell)
resolution : 460, 1090 (x, y)
extent : -3026365, -1546545, 922441.7, 2362332 (xmin, xmax, ymin, ymax)
crs : +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
Sea ice concentration problem:
I wish to use sea ice concentration data from http://cersat.ifremer.fr/data/products/catalogue. Data are available here: (ftp server) or here (3 files I downloaded; < 1 MB each)
This product contains Antarctic polar sterographic 12.5 km resolution grids of sea ice concentration. Every day of the year has its own .nc file. I load 3 files (3 days of sea ice concentration data) in R as follows:
library(raster)
library(ncdf4)
setwd("Computer/covariates/seaice_IFREMER")
#first import all files in a single folder as a list
rastlist <- list.files(path = getwd(), pattern='.nc$', all.files=TRUE, full.names=FALSE)
# stack the raster data
icerasters <- stack(rastlist)
Here I get warnings that
[1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named nj BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
#Plot sea ice data from a single day
plot(icerasters[[1]])
This looks strange! It is a mirror image of Antarctica!
# Use raster::flip to change it:
icerasters <- flip(icerasters, direction='y')
plot(icerasters[[1]])
Better. This looks like the Antarctic Polar Sterographic Projection.
But I still have a problem. The 'values' of the x/y (lat/lon) axis (in other words the extent of the raster) do not correspond to values I associate with the polar stereographic projection.
The units of the polar stereographic projection are meters (just as the units for WGS84 are in degrees). This means the projected bounds (extent) of the polar stereographic projection should be in 'millions', like those of bathymetry in the preamble. This is not what appears on the ice map.
Let us look at the raster information: print(icerasters)
class : RasterStack
dimensions : 664, 632, 419648, 106 (nrow, ncol, ncell, nlayers)
resolution : 1, 1 (x, y)
extent : 0.5, 632.5, 0.5, 664.5 (xmin, xmax, ymin, ymax)
crs : NA
I can set crs manually, since it is not recorded in the imported raster:
Polar.prj = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs "
crs(icerasters) = Polar.prj
But this does not change the values of the extent / axes. I need these to correspond to those mapped in the bathymetry layer.
How do I do that?
The IFREMER sea ice data are derived from NSIDC 12.5 km data. I have exactly the same problem downloading data from NSIDC (https://nsidc.org/data/ae_si12). These data are stored as .he5 files. Whether .nc or .he5, I always get the warning:
**** WARNING **** I was asked to get a varid for dimension named nj BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"
when importing the data from this sea ice product. Perhaps this is the root of the problem.
