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:
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:
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  "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[])
This looks strange! It is a mirror image of Antarctica!
# Use raster::flip to change it: icerasters <- flip(icerasters, direction='y') plot(icerasters[])
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.