Script for generating .kml files for Ocean and Sea Ice Satellite Application Facility (OSI-SAF) products.

OSISAF_kml_generator.py

""" Script for generating .kml files for Ocean and Sea Ice Satellite Application Facility(OSI-SAF) products.

    Made as a show case.


    USAGE:
    python OSISAF_kml_generator.py [with arguments]
    
    EXAMPLES:
    python OSISAF_kml_generator.py -p conc
    python OSISAF_kml_generator.py -p edge -y 2015 -m 04 -d 04 -hem sh -cmap rainbow
    python OSISAF_kml_generator.py -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
    COMMENTS:
    A show case script only supporting three products at the moment.
===============================================
Name:          OSISAF_kml_generator.py
Author(s):     Trygve Halsne 08.02.2018 (dd.mm.YYYY)
Modifications:
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 = datetime.datetime.now()

# 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",
            "type":"sea_ice_classification"}
ODAP_band_list = {"conc":"ice_conc", "edge":"ice_edge",
            "type":"ice_type"}

# Set TDS path
cat = TDSCatalog('http://thredds.met.no/thredds/catalog/osisaf/met.no/ice/{}/{}/{}/catalog.xml'.format(product, year, month))

# Dataset naming convention
myDataSet = str("ice_" + product + '_' + hemisphere + '_'+ projection + '-100_multi_' + year+month+day+'1200.nc')
dataset = cat.datasets[myDataSet]

# Parse OPeNDAP stream
ds = Dataset(dataset.access_urls['OPENDAP'])
max_value=ds.variables[ODAP_band_list[product]][:].max()
min_value=ds.variables[ODAP_band_list[product]][:].min()
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
    times.append(t0+datetime.timedelta(0,t))


# Parse bounding box from WMS
request = "?service=WMS&version=1.3.0&request=GetCapabilities" # URL GetCapabilities
wms = dataset.access_urls['WMS']
f=ul2.urlopen(wms+request)
d = f.read()
f.close()
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]
    nsmap.pop(None)
bb=root.findall('.//default:EX_GeographicBoundingBox',namespaces=nsmap)

# 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'
kml.save(outname)

if os.path.isfile(outname):
    print "Success!\nScript execution time: %s. \nThe file: %s, is written to this folder." % (datetime.datetime.now()-s_t0, outname)
else:
    print "Something happened during script execution. Please debug and find out what caused the problem."