A script for extracting Bedmap2 masks as geometric features.

extractBedmap2Antarctica.py

from netCDF4 import Dataset
import numpy
from matplotlib import pyplot

from shapely.geometry import Polygon, mapping
from shapely.ops import cascaded_union

from utils.feature_write_utils import write_all_features

def extract_geometry(mask):
    mask = numpy.logical_or(mask, bathymetry > 0.)
    expanded_lon = numpy.zeros((lon.size+2))
    expanded_lat = numpy.zeros((lat.size+2))
    expanded_mask = numpy.zeros((mask.shape[0]+2,mask.shape[1]+2))
    expanded_lon[1:-1] = lon
    expanded_lon[0] = 2*lon[0]-lon[1]
    expanded_lon[-1] = 2*lon[-1]-lon[-2]
    expanded_lat[1:-1] = lat
    expanded_lat[0] = 2*lat[0]-lat[1]
    expanded_lat[-1] = 2*lat[-1]-lat[-2]
    expanded_mask[1:-1,1:-1] = mask
    
    
    cs = pyplot.contour(expanded_lon,expanded_lat,expanded_mask,[0.99999])
    pyplot.draw()
    paths = cs.collections[0].get_paths()
    polys = []
    for path in paths:
        v = path.vertices
        x = v[:,0]
        y = v[:,1]
        poly = Polygon([(i[0], i[1]) for i in zip(x,y)])
        polys.append(poly)
        
    return mapping(cascaded_union(polys))

inFile = Dataset('/home/xylar/data/observations/Antarctica/bedmap2_on_ETOPO1_grid.nc')

lon = inFile.variables['lon'][:]
lat = inFile.variables['lat'][:]
print lat
exit()
bathymetry = inFile.variables['bathymetry'][:,:]
ice_mask = inFile.variables['ice_mask'][:,:]
grounded_mask = inFile.variables['grounded_mask'][:,:]

masks = {}
masks['AntarcticIceCoverage'] = ice_mask
masks['AntarcticGroundedIceCoverage'] = grounded_mask
features_file = {}
features_file['features'] = []
for name in masks:
    
    properties = {}
    properties['name'] = name
    properties['component'] = 'bedmap2'
    properties['author'] = 'https://www.bas.ac.uk/project/bedmap-2/'
    properties['object'] = 'region'
    properties['tags'] = ''
    feature = {}
    feature['properties'] = properties
    feature['geometry'] = extract_geometry(masks[name])
    features_file['features'].append(feature)

out_file_name = "features.geojson"

out_file = open(out_file_name, 'w')

out_file.write('{"type": "FeatureCollection",\n')
out_file.write(' "groupName": "enterNameHere",\n')
out_file.write(' "features":\n')
out_file.write('\t[\n')
write_all_features(features_file, out_file, '\t\t')
out_file.write('\n')
out_file.write('\t]\n')
out_file.write('}\n')