Merge branch xylar/add_ice_shelf_feature_script to master

Add scripts for creating ice-shelf features

feature_creation_scripts/README.md

@@ -0,0 +1,3 @@
+This folder is for scripts that were used to create geometric features elsewhere
+in the repository.  They should all contain instructions for how to obtain or generate
+the datasets they use.

feature_creation_scripts/iceshelves/README.md

@@ -0,0 +1,18 @@
+These scripts are used to create the `iceshelf` features that cover not only
+the current positions of ice shelves but also extend them both under grounded
+ice and out to the -800 m isobath, allowing floating regions to be identified
+with ice shelves even after the grounding line has advanced or retreated.
+
+Instructions for re-creating the `iceshelf` feature sin `iceshelves.geojson`:
+
+ * Download `bedmap2_bin.zip` from [https://www.bas.ac.uk/project/bedmap-2/](https://www.bas.ac.uk/project/bedmap-2/)
+ * Unzip `bedmap2_bin.zip` (should create `./bedmpa2_bin/`)
+ * Download `antarctica_ice_velocity_900m.nc` from
+   [http://nsidc.org/data/docs/measures/nsidc0484_rignot/](http://nsidc.org/data/docs/measures/nsidc0484_rignot/)
+ * Do the following:
+```bash
+conda install numpy scipy netcdf4 matplotlib progressbar scikit-fmm basemap \
+    scikit-image shapely descartes cartopy
+ln -s feature_creation_scripts/iceshelves/driver.bash
+./driver.bash
+```

feature_creation_scripts/iceshelves/bedmap2_to_netcdf.py

@@ -0,0 +1,116 @@
+#!/usr/bin/env python2
+
+from netCDF4 import Dataset
+import numpy
+
+from progressbar import ProgressBar, Percentage, Bar, ETA
+"""
+This script converts the Bedmap2 binary dataset into a NetCDF file, including
+longitude and latitude fields.  The Bedmap2 bindary dataset is available via
+https://www.bas.ac.uk/project/bedmap-2/ and directly downloadable from
+https://secure.antarctica.ac.uk/data/bedmap2/bedmap2_bin.zip
+
+The Bedamp2 binary data is assumed to reside in the folder ./bedmap2_bin/.
+This script produces the file bedmap2.nc in the current working directory.
+"""
+
+def readBedmap2Var(name, fillValue = -9999):
+  field = numpy.fromfile('%s/bedmap2_%s.flt'%(inFolder,name),
+               dtype = numpy.float32, count=nx*ny)
+  field = field.reshape(nx,ny)
+  field = field[::-1,:]
+  if fillValue is not None:
+      field = numpy.ma.masked_array(field, field == fillValue)
+
+  return field
+
+def addVariable(field, varName, outType):
+    var = ncfile.createVariable(varName, outType, ('nx','ny'))
+    var[:,:] = field
+    print varName, numpy.amin(field), numpy.amax(field)
+
+
+def makeBedmap2LonLat():
+    projectionLat = -71.0*numpy.pi/180.0
+    a = 6378137.0 # m  equatorial radius of the WGS84 ellipsoid according to Wikipedia
+    f = 1/298.257223563 # oblateness of the WGS84 ellipsoid
+    latToParametricLatFactor = numpy.sqrt(1.0 - (2*f - f**2)) # sqrt(1 - e^2)
+    sinProjectionLat = numpy.sin(numpy.arctan(latToParametricLatFactor*numpy.tan(projectionLat)))
+
+    # find the lon/lat of each bedmap2 point
+    x = dx*(numpy.arange(nx)-0.5*(nx-1))
+    y = dx*(numpy.arange(ny)-0.5*(ny-1))
+    (X, Y) = numpy.meshgrid(x,y,indexing='ij')
+    R = numpy.sqrt(X**2 + Y**2)
+    Lon = numpy.arctan2(Y,X)*180./numpy.pi
+
+    del X, Y
+
+    stretch = 1.0
+
+    iterCount = 30
+    widgets = ['Computing latitudes: ', Percentage(), ' ', Bar(), ' ', ETA()]
+    iterBar = ProgressBar(widgets=widgets, maxval=iterCount ).start()
+
+    for iterIndex in range(iterCount):
+        cosLat = stretch*R/a
+        sinLat = -numpy.sqrt(1 - cosLat**2)
+        stretch = (1.0 - sinLat)/(1.0 - sinProjectionLat)
+        print "sinLat bounds:", numpy.amax(sinLat), numpy.amin(sinLat)
+        print "stretch bounds:", numpy.amax(stretch), numpy.amin(stretch)
+        iterBar.update(iterIndex+1)
+    iterBar.finish()
+
+    mask = cosLat != 0.0
+    tanLat = sinLat[mask]/cosLat[mask]/latToParametricLatFactor
+    Lat = -90.*numpy.ones(Lon.shape)
+    Lat[mask] = numpy.arctan(tanLat)*180./numpy.pi
+
+    return (x, y, Lon, Lat)
+
+inFolder = './bedmap2_bin'
+dx = 1e3 #m
+nx = 6667
+ny = 6667
+
+ncfile = Dataset('Bedmap2.nc', 'w', format='NETCDF4')
+
+ncfile.createDimension('nx',nx)
+ncfile.createDimension('ny',ny)
+
+iceThickness = readBedmap2Var('thickness')
+surface = readBedmap2Var('surface')
+draft = surface - iceThickness
+
+addVariable(surface,'iceSurface','f4')
+addVariable(iceThickness,'iceThickness','f4')
+addVariable(draft,'iceDraft','f4')
+
+del iceThickness, surface, draft
+
+addVariable(readBedmap2Var('bed'),'bed','f4')
+
+bedmap2Mask = readBedmap2Var('icemask_grounded_and_shelves', fillValue = None)
+
+addVariable(bedmap2Mask,'bedmap2Mask','i4')
+
+addVariable(numpy.array(bedmap2Mask < 0., int),'openOceanMask','i4')
+addVariable(numpy.array(bedmap2Mask == 1., int),'floatingIceMask','i4')
+bareBedrockMask = numpy.logical_not(readBedmap2Var('rockmask').mask)
+groundedMask = numpy.logical_and(bedmap2Mask == 0., numpy.logical_not(bareBedrockMask))
+addVariable(numpy.array(groundedMask,int),'groundedIceMask','i4')
+addVariable(numpy.array(bareBedrockMask,int),'bareBedrockMask','i4')
+
+(x, y, Lon, Lat) = makeBedmap2LonLat()
+
+var = ncfile.createVariable('x', 'f4', ('nx',))
+var[:] = x
+var = ncfile.createVariable('y', 'f4', ('ny',))
+var[:] = y
+addVariable(Lon,'lon','f4')
+addVariable(Lat,'lat','f4')
+
+
+ncfile.close()
+
+

feature_creation_scripts/iceshelves/driver.bash

@@ -0,0 +1,50 @@
+#!/bin/bash
+
+# This is the main driver script.  See README.md for details.
+
+rm -f iceshelves.geojson
+
+# convert Bedmap2 data to NetCDF format
+./feature_creation_scripts/iceshelves/bedmap2_to_netcdf.py
+
+# Compute the distance from each point to the closest ice shelf
+./feature_creation_scripts/iceshelves/write_iceshelf_distance.py
+
+# Regrid the Measures velocity data to the Bedmap2 grid
+./feature_creation_scripts/iceshelves/ice_vel_to_bedmap2.py
+
+# Plot streamlines for separating basins on ice shelves
+./feature_creation_scripts/iceshelves/plot_streamline.py
+
+# Write out a map of ice-shelf numbers
+ln -sfn feature_creation_scripts/iceshelves/ice_shelf_lat_lon.csv
+./feature_creation_scripts/iceshelves/write_iceshelf_numbers.py
+
+# Create combine the IMBIE basin features
+./driver_scripts/setup_antarctic_basins.py
+
+# Make images for each IMBIE basin
+./feature_creation_scripts/iceshelves/write_basin_images.py
+
+# Find the distance to each basin
+./feature_creation_scripts/iceshelves/write_basin_distances.py
+
+# Write out a map of basin numbers
+./feature_creation_scripts/iceshelves/write_basin_numbers.py
+
+# Find which shelves overlap with which IMBIE basins
+./feature_creation_scripts/iceshelves/write_expanded_overlaps.py
+
+# Create features for each ice shelf in each basin
+./feature_creation_scripts/iceshelves/write_iceshelf_features.py
+
+# Fix features that cross the antemeridian (+/- 180)
+./fix_features_at_antimeridian.py -f iceshelves_before_fix.geojson \
+    -o iceshelves.geojson
+rm -f iceshelves_before_fix.geojson
+
+# Plot the features
+./plot_features.py -f iceshelves.geojson -m southpole
+
+# To split the features, uncomment the following
+# ./split_features -f iceShelfes.geojson

feature_creation_scripts/iceshelves/ice_shelf_lat_lon.csv

@@ -0,0 +1,101 @@
+Ice Shelf Name,Lat (o),Lon (o),Area
+Larsen G,-74.6,-61.76,412.07
+Larsen F,-74.07,-61.23,827.833
+Larsen E,-73.4,-61.08,1184.3891
+Larsen D,-70.7,-61.72,22548.08
+Larsen C,-67.46,-62.73,46465.49
+Larsen B,-65.55,-61.42,6754.5
+Wordie,-69.25,-67.3,276.85
+Wilkins,-70.49,-72.07,12866.06
+Bach,-72.04,-72.1,4578.5
+George VI,-72.35,-69.62,23433.8
+Stange,-73.22,-76.67,8027.2
+Ronne,-78.45,-63.51,338886.5
+Ferrigno,-73.63,-83.57,117.4
+Venable,-73.09,-87.19,3193.8
+Abbot,-72.79,-96.37,29688
+Cosgrove,-73.53,-100.46,3033.4
+Pine Island,-74.93,-100.72,6249.4
+Thwaites,-75.07,-106.41,5499.02
+Crosson,-75,-110.37,3229.2
+Dotson,-74.62,-112.75,5802.5
+Getz,-74.41,-124.38,34018.0712
+Land,-75.59,-141.41,639.983
+Nickerson,-75.8,-145.81,6495.3
+Sulzberger,-77.09,-148.9,12333.2
+Swinburne,-77.35,-153,899.866
+Withrow,-77.17,-157.08,631.779
+Ross West,-80.82,-167.92,306104.618
+Ross East,-80.48,171.89,194703.9
+Drygalski,-75.39,163.37,2337.956
+Nansen,-74.85,163.25,1985.08
+Aviator,-73.91,165.44,785.2
+Mariner,-73.32,168.08,2705.41
+Lillie,-70.79,164,769.711
+Rennick,-70.73,161.78,3272.9
+Cook,-68.56,152.74,3461.5
+Ninnis,-68.26,147.54,1899.4
+Mertz,-67.26,145.35,5522.4
+Dibble,-66.24,134.74,1481.55
+Holmes,-66.78,127.22,1921.2
+Moscow University,-67.01,120.27,5797.6
+Totten,-67.06,115.7,6032.32
+Vincennes,-66.72,109.95,934.96
+Conger/Glenzer,-65.78,103.41,1547.1
+Tracy/Tremenchus,-65.86,101.21,2844.7
+Shackleton,-65.98,97.52,26079.656
+West,-66.95,85,15665.56
+Publications,-69.76,75.16,1551.27
+Amery,-70.38,70.43,60654
+Wilma/Robert/Downer,-67.03,56.6,857.53
+Edward VIII,-66.73,56.43,411.3
+Rayner/Thyer,-67.73,48.45,640.737
+Shirase,-69.95,38.55,820.62
+Prince Harald,-69.36,36.19,5391.55
+Baudouin,-70.06,28.89,32951.8
+Borchgrevink,-70.26,19.81,21579.8
+Lazarev,-69.93,14.47,8519.48
+Nivl,-70.29,11.11,7284.9557
+Vigrid,-70.24,8.38,2088.6
+Fimbul,-70.36,-0.43,40843.4
+Jelbart,-70.88,-4.66,10843.8
+Atka,-70.63,-6.88,1969.1
+Ekstrom,-71.04,-8.62,6872
+Quar,-71.23,-10.91,2155.6
+Riiser-Larsen,-73.13,-16.48,43450.01
+Brunt/Stancomb,-74.93,-22.61,36893.98
+Filchner,-80.28,-42.29,104253.3
+Mendelssohn,-71.3,-72.76,471.337
+Brahms,-71.47,-73.64,238.161
+Verdi,-71.6,-74.49,110.601
+Hull,-75.09,-136.89,210.824
+Richter,-77.09,204.58,99.4619
+Hamilton,-77.51,202.09,183.672
+Moubray,-71.95,170.4,151.677
+Tucker,-72.64,169.7,443.365
+Dennistoun,-71.18,168.05,88.0504
+SmithInlet,-70.99,167.58,92.9502
+Wylde,-73.62,166.84,124.279
+Fitzgerald,-73.68,166.44,316.159
+Tinker,-74.1,165.02,110.538
+Campbell,-74.62,164.42,71.1901
+HarbordGlacier,-75.93,162.76,62.8208
+Cheetham,-75.75,162.73,98.8261
+GeikieInlet,-75.58,162.63,273.471
+Nordenskjold,-76.2,162.58,236.545
+Suvorov,-69.96,160.4,216.277
+Gillet,-69.56,159.71,123.939
+Noll,-69.41,159.07,108.881
+Matusevitch,-69.22,157.36,517.311
+Slava,-68.75,154.87,275.682
+PourquoiPas,-66.23,135.72,188.818
+Frost,-67.07,128.54,172.08
+HolmesWest,-66.36,126.24,340.072
+Voyeykov,-66.63,124.52,671.106
+Underwood,-66.7,107.88,161.995
+Helen,-66.64,93.88,290.611
+Utsikkar,-67.54,61.25,92.2136
+Mulebreen,-67.47,59.35,297.108
+Myers,-67.21,49.81,103.162
+Zubchatyy,-67.24,49.03,269.43
+Hannan,-67.52,47.39,212.079

feature_creation_scripts/iceshelves/ice_vel_to_bedmap2.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env python
+
+import numpy
+from netCDF4 import Dataset
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+
+from mpl_toolkits.basemap import interp
+
+import scipy.ndimage.filters as filters
+
+inFileName = 'Bedmap2_ice_shelf_distance.nc'
+inFile = Dataset(inFileName,'r')
+X = inFile.variables['X'][:,:]
+Y = inFile.variables['Y'][:,:]
+inFile.close()
+
+velocityFileName = 'antarctica_ice_velocity_900m.nc'
+inFile = Dataset(velocityFileName, 'r')
+vx = inFile.variables['vx'][::-1,:]
+vy = inFile.variables['vy'][::-1,:]
+inFile.close()
+
+
+(ny,nx) = vx.shape
+
+dx = 9e2
+xMin = -2800000.0
+yMax = 2800000.0
+
+x = xMin + dx*numpy.arange(nx)
+y = yMax + dx*(numpy.arange(ny)-(ny-1))
+
+maskValue=-1e34
+Vx = interp(vx,x,y,X,Y,masked=maskValue)
+Vx[Vx == maskValue] = 0.0
+Vy = interp(vy,x,y,X,Y,masked=maskValue)
+Vy[Vy == maskValue] = 0.0
+
+mask = 1.0*numpy.logical_or(vx != 0, vy != 0)
+Mask = interp(mask,x,y,X,Y,masked=maskValue)
+Mask[Mask == maskValue] = 0.0
+
+(ny,nx) = Vx.shape
+
+filterSigma = 10.0
+
+Vx = filters.gaussian_filter(Vx,filterSigma)
+Vy = filters.gaussian_filter(Vy,filterSigma)
+Mask = filters.gaussian_filter(Mask,filterSigma)
+
+mask = Mask > 1e-6
+Vx[mask] /= Mask[mask]
+Vy[mask] /= Mask[mask]
+
+
+outFileName = 'Bedmap2_grid_velocities.nc'
+outFile = Dataset(outFileName,'w')
+
+outFile.createDimension('x',nx)
+outFile.createDimension('y',ny)
+var = outFile.createVariable('vx','f8',('y','x'))
+var[:,:] = Vx
+var = outFile.createVariable('vy','f8',('y','x'))
+var[:,:] = Vy
+
+outFile.close()
+
+#extent = [xMin,numpy.amax(x),yMax,numpy.amin(y)]
+
+#extent = [numpy.amin(X),numpy.amax(X),numpy.amax(Y),numpy.amin(Y)]
+#
+#plt.figure()
+#plt.imshow(Vx,extent=extent,interpolation='nearest')
+#plt.colorbar()
+#plt.gca().invert_yaxis()
+#plt.figure()
+#plt.imshow(Vy,extent=extent,interpolation='nearest')
+#plt.colorbar()
+#plt.gca().invert_yaxis()
+#plt.figure()
+#plt.imshow(Mask,extent=extent,interpolation='nearest')
+#plt.colorbar()
+#plt.gca().invert_yaxis()
+#plt.show()