Sea Ice Frequency - npolar/RemoteSensing GitHub Wiki

Explanation of the script, which is found here

OUT OF DATE .... NOT UPDATED YET

Short description of the steps:

  1. Import all NSIDC charts into GeoTIFF and map project

    (At present daily images taken as is, no 5-day median yet)

  2. Ice Persistence Map:

    For every file in the requested time:

    If ice is more than 15%: add (100 / NumberOfDays) in ice_persistency raster

    Otherwise: leave as is

  3. Min/Max map:

    a) Create map showing total number of ice-days

    b) Maximum map: where ice-day = 0 (never ever ice) set 0, else set 1 (at least one day of ice)

    c) Minimum map: where ice-days = NumberOfDays (absolute every day ice) set 1, otherwise 0

  4. Clean up noise: Sieve Filter removing singular pixels (less then 2 neighbouring ones=

  5. Correct errors at coast:

    a) VARIANT ONE: In problematic areas only, if ice less than 1.5% of total days (30 years = less than 13 days), remove these. (Extrem winters with long duration will be caught)

    b) VARIANT TWO: Count as ice if 2 days before and after also ice present. (Not Programmed Yet).

Core of Program

Create filelist including all files for the given month between startyear and stopyear inclusive, then the the following functions follow (see code for details):

def Bin2GeoTiff(infile,outfilepath )

This function takes the NSIDC charts, being a flat binary string, and converts them to GeoTiff. Some details here: Reading NSIDC Files

Info on the ice concentration charts: http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
Info on the map projection: http://nsidc.org/data/polar_stereo/ps_grids.html

The GeoTiff files are map projected to EPSG:3411, being the NSIDC-specific projection. There also exists a GeoTiff reprojected to EPSG:3575 which is the NP-standard for Barents/Fram-Strait. Details on how to map project are found here: Map Projection NSIDC charts Map Projected NSIDC chart Projected to EPSG3575

def CreateIcePercistanceMap(inpath, outfilepath)

This function creates the ice persistence charts from the sea ice concentration charts. The function loops through each concentration map, if the value is larger than 38 = 15.2%, the value of "100.0 / NumberOfDays" is added --> if there is ice every single day, that pixel will be "100"

outarray = numpy.where( (iceraster >=  38), (outarray + ( 100.0 / NumberOfDays ) ) , outarray)
outarray = numpy.where( (iceraster ==   251), 251 , outarray)
outarray = numpy.where( (iceraster ==   252), 252 , outarray)
outarray = numpy.where( (iceraster ==   253), 253 , outarray)
outarray = numpy.where( (iceraster ==   254), 254 , outarray)
outarray = numpy.where( (iceraster ==   255), 255 , outarray)

The output is also available as EPSG:3575

def CreateMaxMinIce(inpath, outfilepath):

This function creates a minimum and maximum map, and as "side-effect" a map showing number of days with ice. The output is available as GeoTIFF and as shapefile, both in EPSG:3411 and EPSG:3575.

 * maximum = at least one day ice at this pixel
 * minimum = every day ice at this pixel
 * In addition a file simply giving the number of days with ice
     
 * The poly shapefile has all features as polygon, the line shapefile only the max or min ice edge
 * The shapefiles select only the largest polygon such that noise along coast-lines is reduced and
   only the main ice edge shown. Check size value, may need to be adjusted for low ice years when 
   largest polygon very small(?).

image

image

def EPSG3411_2_EPSG3575(infile)

a function reprojecting from EPSG:3411 to EPSG:3575

def ReprojectShapefile(infile, inproj = "EPSG:3411", outproj = "EPSG:3575")

a function reprojecting shapefiles from EPSG:3411 to EPSG:3575

Removing noise and erroneous areas

THe figure above shows some singular pixels around the fast ice in the minimum map. This is takenm away by the gdal Sieve filter:

gdal.SieveFilter( srcband, maskband, dstband,threshold = 3, connectedness = 4 )

More difficult is the misclassification along the coast lines of Norway and Iceland and others. The suggested solution uses a mask, but examines if ice was extended periods within the mask to account for real ice, especially at the edge of the mask:

Step 1: This masks indicates the problematic areas -- it is buffered along the "problem coasts":

Step 2: Pixels within this mask, i.e. along the problematic coasts, are only removed when less ice cover is less than 1.5% of the total days (10 yer period = 300 total days, 20 year period = 900 total days).

Step 3: A sieve filter filtering within this mask singular pixels. For March 1980-2010 the yellow area is removed in this process:

For March 1990-2010 the pink area is thus removed:

Questions:

Should min/max be calculated fro EPSG3575 from reprojected raster or from original one?
Smooth ice edge shape files to not show single pixels?
Wrong areas -- remove polygons where DN = 1 and size lesser then x

###comments add area for shapes http://gis.stackexchange.com/questions/83869/is-there-any-way-to-add-an-area-column-to-new-shapefile-created-from-raster-with
smooth http://stackoverflow.com/questions/1849928/how-to-intelligently-degrade-or-smooth-gis-data-simplifying-polygons

-- thoughts on removing noise: Apply gdal sieve and then(!) burn in landmask since sieve changes coastline -- for maximum chart -- reraster polygon

⚠️ **GitHub.com Fallback** ⚠️