icebed.R

file <- "~/Git/icebed/BedMachineGreenland-2017-09-20.nc"
## to GDAL a NetCDF variable is a SUBDATASET
# Subdatasets:
# SUBDATASET_1_NAME=NETCDF:"BedMachineGreenland-2017-09-20.nc":mask
# SUBDATASET_1_DESC=[18346x10218] mask (8-bit integer)
# SUBDATASET_2_NAME=NETCDF:"BedMachineGreenland-2017-09-20.nc":surface
# SUBDATASET_2_DESC=[18346x10218] surface_altitude (16-bit integer)
# SUBDATASET_3_NAME=NETCDF:"BedMachineGreenland-2017-09-20.nc":thickness
# SUBDATASET_3_DESC=[18346x10218] land_ice_thickness (16-bit integer)
# SUBDATASET_4_NAME=NETCDF:"BedMachineGreenland-2017-09-20.nc":bed
# SUBDATASET_4_DESC=[18346x10218] bedrock_altitude (16-bit integer)
# SUBDATASET_5_NAME=NETCDF:"BedMachineGreenland-2017-09-20.nc":errbed
# SUBDATASET_5_DESC=[18346x10218] errbed (16-bit integer)
# SUBDATASET_6_NAME=NETCDF:"BedMachineGreenland-2017-09-20.nc":source
# SUBDATASET_6_DESC=[18346x10218] source (8-bit integer)
# SUBDATASET_7_NAME=NETCDF:"BedMachineGreenland-2017-09-20.nc":geoid
# SUBDATASET_7_DESC=[18346x10218] geoid_height_above_reference_ellipsoid (16-bit integer)
# 

library(raster)  ## drives NetCDF via the R package ncdf4

## Loading required package: sp

## to raster, a 2D slice from a variable is a band analogous to a GDAL 2D dataset
surface <- raster(file, varname = "surface")

## Loading required namespace: ncdf4

## Warning in .getCRSfromGridMap4(atts): cannot process these parts of the CRS:
## geoid=eigen-6c4

## raster gets the projection wrong, as it misses the standard parallel
projection(surface) <- "+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=-45 +lat_ts=70 +ellps=WGS84"

## extent in projection
extent(surface)

## class       : Extent 
## xmin        : -653000 
## xmax        : 879700 
## ymin        : -3384500 
## ymax        : -632600

projection(surface)

## [1] "+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=-45 +lat_ts=70 +ellps=WGS84"

## get a map
library(rnaturalearth)
library(rnaturalearthhires)
library(sf)

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

coast <- st_transform(ne_coastline(scale = 10, returnclass = "sf"), projection(surface))
plot(surface, col = viridis::viridis(64))
plot(st_geometry(coast), add = TRUE)

## but GDAL gets it right, though I need to hack the sds name ...
gdal <- rgdal::readGDAL('NETCDF:"BedMachineGreenland-2017-09-20.nc":surface')

## NETCDF:"BedMachineGreenland-2017-09-20.nc":surface has GDAL driver netCDF 
## and has 18346 rows and 10218 columns

sp::proj4string(gdal)

## [1] "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

sp::bbox(gdal)

##        min     max
## x  -653000  879700
## y -3384500 -632600