Skip to content

Instantly share code, notes, and snippets.

Last active February 8, 2018 12:35
Show Gist options
  • Save hevgyrt/b507d0405f62086770f083bb8957e747 to your computer and use it in GitHub Desktop.
Save hevgyrt/b507d0405f62086770f083bb8957e747 to your computer and use it in GitHub Desktop.
Script for generating .kml files for Ocean and Sea Ice Satellite Application Facility (OSI-SAF) products.
""" Script for generating .kml files for Ocean and Sea Ice Satellite Application Facility(OSI-SAF) products.
Made as a show case.
python [with arguments]
python -p conc
python -p edge -y 2015 -m 04 -d 04 -hem sh -cmap rainbow
python -p edge -y 2016 -m 12 -d 14 -hem nh
required arguments:
-p, --product PRODUCT
Product Type. Supporting "edge", "conc", "type"
optional arguments:
-h, --help show this help message and exit
-y, --year YEAR Year, e.g. "2017"
-m, --month MONTH
Month, e.g. "03" (March)
-d, --day DAY
Day, e.g. "01"
-hem, --hemisphere HEMISPHERE
Hemisphere ("nh", "sh")
-cmap, --colormap COLORMAP
Colormap for output. Defualt greyscale. See
GetCapabilities for more options
A show case script only supporting three products at the moment.
Author(s): Trygve Halsne 08.02.2018 (
Copyright: (c) Norwegian Meteorological Institute, 2018
import datetime
import simplekml
from netCDF4 import Dataset
import numpy as np
from siphon.catalog import TDSCatalog
import urllib2 as ul2
from lxml import etree
import argparse
import os
s_t0 =
# Read input values from command line using argparser
parser = argparse.ArgumentParser(description='Generate .kml file for OSI-SAF ice-edce, ice-concentration and ice-type products.')
parser.add_argument("-p", "--product", help='Product Type. Supporting "edge", "conc", "type"', type=str, required=True)
parser.add_argument("-y", '--year', help=('Year, e.g. "2017"'), type=str , default="2017")
parser.add_argument("-m", '--month', help=('Month, e.g. "03" (March)'),type=str , default="01")
parser.add_argument("-d", '--day', help=('Day, e.g. "01"'),type=str , default="13")
parser.add_argument("-hem", '--hemisphere', help=('Hemisphere ("nh", "sh")') ,type=str , default="nh")
parser.add_argument("-cmap", '--colormap', help=('Colormap for output. Defualt greyscale. See GetCapabilities for more options'),type=str , default="greyscale")
args = parser.parse_args()
product = getattr(args,'product')
year = getattr(args,'year')
month = getattr(args,'month')
day = getattr(args,'day')
hemisphere = getattr(args,'hemisphere')
colormap = getattr(args,'colormap')
# Set some initial variables
projection = "polstere" #ease?
dateformat = "%Y-%m-%dT%H:%M:%S.%fZ"
band_list = {"conc":"sea_ice_area_fraction", "edge":"sea_ice_classification",
ODAP_band_list = {"conc":"ice_conc", "edge":"ice_edge",
# Set TDS path
cat = TDSCatalog('{}/{}/{}/catalog.xml'.format(product, year, month))
# Dataset naming convention
myDataSet = str("ice_" + product + '_' + hemisphere + '_'+ projection + '-100_multi_' + year+month+day+'')
dataset = cat.datasets[myDataSet]
# Parse OPeNDAP stream
ds = Dataset(dataset.access_urls['OPENDAP'])
time = ds.variables['time'][:]
lonmin = ds.variables['lon'][:].min()
lonmax = ds.variables['lon'][:].max()
latmin = ds.variables['lat'][:].min()
latmax = ds.variables['lat'][:].max()
t0=datetime.datetime(1978,1,01, 00,00,00)
times = []
for t in time: # in case of temporal products
# Parse bounding box from WMS
request = "?service=WMS&version=1.3.0&request=GetCapabilities" # URL GetCapabilities
wms = dataset.access_urls['WMS']
d =
root = etree.fromstring(d,base_url=str(wms+request))
nsmap = root.nsmap
# Tweak the nsmap in order to use ET find and findall
if nsmap.has_key(None):
nsmap['default'] = nsmap[None]
# create the .kml file
kml = simplekml.Kml()
doc = kml.newfolder(name = band_list[product])
for i in range(len(times)):
overlay = doc.newgroundoverlay(name=ODAP_band_list[product])
tmp_time = datetime.datetime.strftime(times[i],dateformat)
href = wms + '?VERSION=1.1.1&REQUEST=GetMap&bbox=%f,%f,%f,%f&SRS=EPSG:4326&WIDTH=2048&HEIGHT=2048&LAYERS=%s&STYLES=boxfill/%s&TRANSPARENT=TRUE&FORMAT=image/png&COLORSCALERANGE=%f,%f&NUMCOLORBANDS=253&TIME=%s' %(lonmin,latmin,lonmax,latmax,ODAP_band_list[product],colormap, min_value,max_value, tmp_time)
overlay.icon.href = href
overlay.latlonbox.west = lonmin
overlay.latlonbox.south = latmin
overlay.latlonbox.east = lonmax
overlay.latlonbox.north = latmax
overlay.icon.refreshmode ='onStop'
overlay.icon.viewboundscale = 0.75
overlay.timestamp = simplekml.TimeStamp(times[i])#.isoformat()+'Z')
#save the kml
outname = ODAP_band_list[product] +'_'+tmp_time.replace(':','-').split('.')[0] + '.kml'
if os.path.isfile(outname):
print "Success!\nScript execution time: %s. \nThe file: %s, is written to this folder." % (, outname)
print "Something happened during script execution. Please debug and find out what caused the problem."
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment