IceChartProcessing - npolar/RemoteSensing GitHub Wiki

== General Info ==

  • ice charts from met.no are downloaded automatically to \henry\isdata (Are maintains the ftp download)
  • the barents_shp folder is the one we are interested in
  • Charts are daily since 1998, before weekly

Ice class text values contain numerical values meaning

  • Open Water = 5%
  • Very Open Drift Ice = 25%
  • Open Drift Ice = 55%
  • Close Drift Ice = 80%
  • Very Close Drift Ice = 95%
  • Fast Ice = 100%

== Limitations of the Data ==

'''Important'''

This data is meant for daily sea ice navigation, taking any source the analyst may find, from 25km passive microwave to 20 meter SAR data, also information a ship may provide. It is not documented which area is based on which maps and it is drawn by visual analysis. It is not a scientific dataset in the sense that results could be reproduced or individual maps reanalysed or strange features in a certain daily map be explained later.

== Software ==

The code is in Python using gdal /ogr libraries (installation inctructions [http://here http://geoinformaticstutorial.blogspot.no/2013/03/installing-python-with-gdal-on-windows.html] and calling the [http://Nest%20Toolbox http://nest.array.ca/web/nest]

It is not using ArcGIS python code in order to be independent of commercial software.

== Processing Steps of Script IceChart.py ==

original code at https://github.com/npolar/RemoteSensing/blob/master/IceChartProcessing.py

Processing Steps:
Reproject Shapefile --> Convert Shape to Raster --> Add Missing Days --> Do Calculations

== met.no Shapefiles ==

The shapefile are unprojected WGS 84 lat/long (EPSG 4326)

Fil:icechart.png

== def ReprojectShapefile(infile) ==

reproject original shapefile to map projection, like EPSG 3575 WGS 84 / North Pole LAEA Europe

open issues:

  • old files have different format, cant be processed yet

Fil:IcechartEPSG3575.png

== def Shape2Raster(shapefile) ==

convert reprojected shapefile to raster -- allows to take all Barentssea or just subset, resolution and extent can be defined.
Uses independently created mask to mask out land areas:

Fil:dailyicechart.jpg

== def AddMissingDays() ==

using the date.time() module:

#http://pymotw.com/2/datetime/
#http://docs.python.org/2/library/datetime.html#strftime-strptime-behavior
d1 =  datetime.date(year1, month1, day1)
d2 =  datetime.date(year2, month2, day2)
diff = datetime.timedelta(days=1)

== def CreateMapFastIceDays () ==

Creates Map where number indicates days with fast ice, 999 is land -- '''not''' considering if days are consecutive

    #Count pixel, if fastice add one, of not keep value -- add land mask
    outarray = numpy.where( (iceraster == 100), outarray + 1 , outarray)
    outarray = numpy.where( (iceraster == 999), 999 , outarray)

Fil:nonconsecutivefastice.jpg

== def CreateMapConsecutiveFastIceDays() ==

Creates Map where number indicates days with fast ice, 999 is land '''considering if days are consecutive'''

Number may be several consecutive periods added or only one -- not differentiated yet

Calculation Method at present. Loop through files and:

    #Read input raster into array
    iceraster = icechart.ReadAsArray()
    
    #Add one where fast ice otherwise keep number, count in countingraster
    outarray = numpy.where( (iceraster == 100), outarray + 1 , outarray)
    countingraster = numpy.where( (iceraster == 100), countingraster + 1 , countingraster)
    
    #Reset pixel to zero when previous day was not fast ice AND if less than wanted number (otherwise pixel is valid)
    outarray = numpy.where( (numpy.logical_and((previousraster == 0) , (countingraster < consecutivenumber))), 0 , outarray)
    #NUMBER IS NOT CONSECUTIVE BUT SEVERAL PERIODS OF CONSECUTIVE ADDED!!
    #reset counting raster to 0 where outarray was set to 0 in line above
    countingraster = numpy.where( (outarray == 0), 0 , countingraster)
    #add landmask
    outarray = numpy.where( (iceraster == 999), 999 , outarray)
           
    #New previousraster
    previousraster = numpy.where( (iceraster == 100), 1 , 0)

after the last round run that because the last pixel has to be compared still:

#Has to finish after last run checking last image
outarray = numpy.where((countingraster < consecutivenumber), 0 , outarray)
outarray = numpy.where( (iceraster == 999), 999 , outarray)

Fil:consecutivefastice.jpg

== def CreatePercentageMap() ==

Creates map showing percentage ice coverage over a given period

Biggest change: Loop replaced by Raster Processing:

#Obsolete pixel-by-pixel version, takes very long

for i in range(rows):

for j in range(cols):

if iceraster[i,j] == 0:

outarray[i,j] = outarray[i,j] + ( 0.0 / NumberOfDays)

elif iceraster[i,j] == 5:

outarray[i,j] = outarray[i,j] + ( 5.0 / NumberOfDays)

elif iceraster[i,j] == 25:

outarray[i,j] = outarray[i,j] + ( 25.0 / NumberOfDays)

elif iceraster[i,j] == 55:

outarray[i,j] = outarray[i,j] + ( 55.0 / NumberOfDays)

elif iceraster[i,j] == 80:

outarray[i,j] = outarray[i,j] + ( 80.0 / NumberOfDays)

elif iceraster[i,j] == 95:

outarray[i,j] = outarray[i,j] + ( 95.0 / NumberOfDays)

elif iceraster[i,j] == 100:

outarray[i,j] = outarray[i,j] + (100.0 / NumberOfDays)

if iceraster[i,j] == 999:

outarray[i,j] = 999.0

      #Array calculation with numpy -- much faster
      outarray = numpy.where( (iceraster ==   0.0), (outarray + (  0.0 / NumberOfDays)) , outarray)
      outarray = numpy.where( (iceraster ==   5.0), (outarray + (  5.0 / NumberOfDays)) , outarray)
      outarray = numpy.where( (iceraster ==  25.0), (outarray + ( 25.0 / NumberOfDays)) , outarray)
      outarray = numpy.where( (iceraster ==  55.0), (outarray + ( 55.0 / NumberOfDays)) , outarray)
      outarray = numpy.where( (iceraster ==  80.0), (outarray + ( 80.0 / NumberOfDays)) , outarray)
      outarray = numpy.where( (iceraster ==  95.0), (outarray + ( 95.0 / NumberOfDays)) , outarray)
      outarray = numpy.where( (iceraster == 100.0), (outarray + (100.0 / NumberOfDays)) , outarray)
    outarray = numpy.where( (iceraster == 999.0), 999.0 , outarray)

Fil:percentage.jpg

== def CreateIceEdgeMap() ==

Create Ice Edge Map

Treshhold the chosen percentage per day.

Example for 30%:

Fil:iceedge30.jpg

== Next Steps: Clip ice edge to get single line ==

manual steps:

  • iceadge --> do polygon to lines

  • the clipping file (iceshape_3575.shp) is polygon

  • Select in iceedge file the iceedge itself

  • Run Geoprocessing Tools>Clip and get the iceedge itself

  • Run again to with selected coast line to remove coast lines from ice edge.

Use this polygon to clip out only lower line of polygon as single line:

Fil:clippinglayer.jpg

Result:

Fil:iceedgelineonly.jpg

== Next Steps: Decadal Average Map ==

A ten year mean is for example taking all February-months 2001-2010 and calculating the average sea ice coverage per pixel using all the ten Februaries.

== Product List ==

BarentsPortal

Wanted by Jon Richard:

Isgrenser:
Siste års isgrense ved 15 % iskonsentrasjon mars og september
Ti års gjennomsnittlig isgrense mars og september
30 års gjennomsnitt isgrense mars og september

Iskonsentrasjonen:
Siste års iskonsentrasjon mars og september*
Ti års gjennomsnittlig iskonsentrasjon mars og september*
30 års gjennomsnitt iskonsentrasjon mars og september*

*Konsentrasjonene er:
80-100% ice concentration 
60-80 
40-60
20-40 
< 20 

Isdekt areal:
Isdekt areal innenfor Barentshavets anerkjente grenser i km2. Dette er en parameter vi vil ta opp med Olga.

Issues:

Barentsportal (John Richard) has ice charts with the classes:

80-100%
60-80%
40-60%
20-40%
<20%

But the met.no data contains according to you info

  • Open Water = 5%
  • Very Open Drift Ice = 25%
  • Open Drift Ice = 55%
  • Close Drift Ice = 80%
  • Very Close Drift Ice = 95%
  • Fast Ice = 100%

So I cannot get the boundaries as in Barentsportal. Could it be that this is a wrong assignment to the met.no-classes (Very open Drift ice being assigned <20% ? ), or how did Harvey get these first thresholds?

== Information from Ola Brandt: ==

Mission: Use the met.no data to create a sea ice indicator index for Svalbard.

Execution: The shape files provided by met.no is rasterized in MATLAB. A binary raster is created comprising 'fast ice' / 'everything else'. The number of days with fast ice is calculated for each pixel and for three different epochs, 1 Sept=>1 Sept (one year), 1 Sept=>1 Mars, 1 Mars=>1 Sept. The latter to indicate inter seasonal differences. Polygons for various regions are applied and the total number of [km^2 days] of ice for each region and epoch are calculated.

Met.no Ice Chart Shape Files are downloaded automatically to \henry\isdata\2012\barents_shp

The following prototype scripts are used:

IceChart2RasterUTM.m: Converts shape files into UTM33 Raster.

IceChartsGiveMOSJProd.m: Calculates number of days with ice for each pixel, filter, splits into various epochs and regions and calculates final [km^2 days]. This script is run in several steps (break command in file).

Follow and read comments within scripts. Note: paths to data and where to store final and intermediate products must be changed.

The MATLAB scripts will call a few sub functions that needs to be accesable:

BatchReaderEndName.m

wgs2utm_fixed_zone.m

As well as data file:

SvalbardBW.mat

Polygon now used are:

C:\Users\ola\Documents\Icecharts\VSpetsPoly.mat'

C:\Users\ola\Documents\Icecharts\NSpetsPoly.mat'

C:\Users\ola\Documents\Icecharts\ESpetsPoly.mat'

C:\Users\ola\Documents\Icecharts\NAustLandPoly.mat'

== Old Obsolete Steps Arcus maps ==

PRELIMINARY AND TOO MANUAL:

  1. Copy the shapefiles from month to C:\Users\max\Documents\IceCharts
  2. If first and last of month missing, consider replacing manually with nearest one
  3. make sure C:\Users\max\Documents\IceCharts\Arcus\EPSG3575 exists and is empty
  4. Run https://github.com/npolar/RemoteSensing/blob/master/IceChart_Arcus.py
  5. Result arcus_map.shp and arcus_map.tif
  6. Repeat for the same month back to 2009
  7. Open arcus_map.shp in QGIS for all years
  8. adjust colours as in, for example, JulyMaps09_13.qgs -- use transparencies and advanced>symbol levels
  9. Use print composer to save as jpg
    10)Open arcus_map.shp from present year and decadal means 30% found at \henry\felles\Forvalt\Miljødata\Harvey\IsAnalyse\Isdata\Mean\Shapefiles
    11)adjust colours as in, for example, DecadalMeans_July.qgs -- use transparencies and advanced>symbol levels
    12)use gimp to crop image -- since areal polygons a not only southern line as extent, artefacts at Northern boundaries

Harveys isanalyse copied for backup to \henry\felles\Forvalt\Miljødata\Max

== Old Working Log ==

Add a custom footer
⚠️ **GitHub.com Fallback** ⚠️