From bdec9f5212e0434ee43a93105e94ea90455e8240 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 25 May 2016 07:35:47 +0000 Subject: [PATCH 01/27] mpop/satout/mitiff.py first draft --- mpop/satout/mitiff.py | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 mpop/satout/mitiff.py diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py new file mode 100644 index 00000000..0d9cf108 --- /dev/null +++ b/mpop/satout/mitiff.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2016. + +# Author(s): + +# Trygve Aspenes + +# This file is part of mpop. + +# mpop is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# mpop is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License along with +# mpop. If not, see . + + +"""mpop mitiff writer interface. +""" + +__revision__ = 0.1 + +import numpy as np +import logging + +logger = logging.getLogger(__name__) + +def save(scene, filename, compression=True, dtype=np.int16, band_axis=2): + """Saves the scene as a MITIFF file. This is local for METNO, but can be read by DIANA + + """ + return mitiff_writer(filename, CFScene(scene, dtype, band_axis), compression=compression) + + +def mitiff__writer(filename, root_object, compression=True): + """ Write data to MITIFF file. """ + From 7f1ca6d6f3b3ba1e9f26ec832d4e02085d69f90e Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Mon, 30 May 2016 09:52:21 +0200 Subject: [PATCH 02/27] Write files ok. Added reading of xml config file, work in progress remove reading of cfg style config --- mpop/satout/mitiff.py | 165 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 160 insertions(+), 5 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 0d9cf108..c3c6c814 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -1,6 +1,8 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # Copyright (c) 2016. +#from tifffile.tifffile import TIFF_TAGS +from TiffImagePlugin import IMAGEDESCRIPTION # Author(s): @@ -26,18 +28,171 @@ __revision__ = 0.1 +import os import numpy as np import logging +from copy import deepcopy +import xml_read +#import xml.etree.ElementTree as etree logger = logging.getLogger(__name__) -def save(scene, filename, compression=True, dtype=np.int16, band_axis=2): - """Saves the scene as a MITIFF file. This is local for METNO, but can be read by DIANA +#------------------------------------------------------------------------------- +# +# Read MiTIFF products config file. +# +#------------------------------------------------------------------------------- +def get_product_config(product_name, force_read=False): + """Read MiTIFF configuration entry for a given product name. + + :Parameters: + product_name : str + Name of MiTIFF product. + + :Arguments: + force_read : Boolean + Force re-reading config file. + + **Notes**: + * It will look for a *mitiff_products.cfg* in MPOP's + configuration directory defined by *PPP_CONFIG_DIR*. + * As an example, see *mitiff_products.cfg.template* in + MPOP's *etc* directory. + """ + return ProductConfigs()(product_name, force_read) + + +class _Singleton(type): + def __init__(cls, name_, bases_, dict_): + super(_Singleton, cls).__init__(name_, bases_, dict_) + cls.instance = None + + def __call__(cls, *args, **kwargs): + if cls.instance is None: + cls.instance = super(_Singleton, cls).__call__(*args, **kwargs) + return cls.instance + +class ProductConfigs(object): + __metaclass__ = _Singleton + + def __init__(self): + self.read_config() + + def __call__(self, product_name, force_read=False): + if force_read: + self.read_config() + return self._products[product_name] + + @property + def product_names(self): + return sorted(self._products.keys()) + + def read_config(self): + from ConfigParser import ConfigParser + + def _eval(val): + try: + return eval(val) + except: + return str(val) + + filename = self._find_a_config_file() + logger.info("Reading MiTIFF config file: '%s'" % filename) + + cfg = ConfigParser() + cfg.read(filename) + products = {} + for sec in cfg.sections(): + prd = {} + for key, val in cfg.items(sec): + prd[key] = _eval(val) + logger.debug("%s %s" % (key,prd[key])) + products[sec] = prd + self._products = products + + @staticmethod + def _find_a_config_file(): + name_ = 'mitiff_products.cfg' + home_ = os.path.dirname(os.path.abspath(__file__)) + penv_ = os.environ.get('PPP_CONFIG_DIR', '') + for fname_ in [os.path.join(x, name_) for x in (home_, penv_)]: + if os.path.isfile(fname_): + return fname_ + raise ValueError("Could not find a MiTIFF config file") + +def save(scene, filename, product_name=None): + """Saves the scene as a MITIFF file. This format can be read by DIANA """ - return mitiff_writer(filename, CFScene(scene, dtype, band_axis), compression=compression) + print "inside save ... " + return mitiff_writer(filename, scene, mitiff_product_name=product_name) -def mitiff__writer(filename, root_object, compression=True): +def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=None): """ Write data to MITIFF file. """ - + print "Inside mitiff_writer ... " + + + name_ = 'mitiff_products.xml' + penv_ = os.environ.get('PPP_CONFIG_DIR', '') + fname_ = os.path.join(penv_,name_) + if os.path.isfile(fname_): + config = xml_read.parse_xml(xml_read.get_root(fname_)) + #raise ValueError("Could not find a MiTIFF config file ", fname_) + + + #if mitiff_product_name: + # options = deepcopy(get_product_config(mitiff_product_name)) + #else: + # print "No MiTIFF product name given. Skip Config. Sure you want this?" + # options = {} + + + #print root_object + tifargs = dict() + args={} + + from libtiff import TIFF + + tif = TIFF.open(filename, mode ='w') + #print dir(tif) + + #np.set_printoptions(threshold=np.nan) + + #tif.TIFFSetField(tif,IMAGEDESCRIPTION,"test") + tif.SetField(IMAGEDESCRIPTION, "TEST") + for ch in root_object.channels: + print ch + print ch.info.get('units', 'None') + #print ch.data + data=ch.data.astype(np.uint8) + tif.write_image(data,) + + tif.close + +if __name__ == '__main__': + from mpop.satellites import PolarFactory + from mpop.projector import get_area_def + import datetime + + logging.basicConfig(level=logging.DEBUG) + + print "Test saving mitiff." + orbit = "37582" + time_slot = datetime.datetime(2016, 5, 24, 15, 25) + global_data = PolarFactory.create_scene("noaa","19","avhrr/3",time_slot, orbit) + global_data.load() + + from pprint import pprint + #pprint(vars(global_data)) + + local_data = global_data.project("eurol_metno", mode="nearest") + + pprint(vars(local_data)) + + save(local_data,"test.mitiff",product_name="AVHRR") + print "Complete." + + + + \ No newline at end of file From c5aa600943b65962d62669e8f66ae8193e58d4c4 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Mon, 30 May 2016 10:05:03 +0200 Subject: [PATCH 03/27] xml_read taken from trollduction and must be adapted --- mpop/satout/xml_read.py | 230 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 230 insertions(+) create mode 100644 mpop/satout/xml_read.py diff --git a/mpop/satout/xml_read.py b/mpop/satout/xml_read.py new file mode 100644 index 00000000..c729b90f --- /dev/null +++ b/mpop/satout/xml_read.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2014, 2015 + +# Author(s): + +# Panu Lahtinen +# Martin Raspaud + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +'''XML reader for Trollduction system and product configuration files. +''' + +import xml.etree.ElementTree as etree +import os +import logging + +LOGGER = logging.getLogger(__name__) + + +class InfoObject(object): + """InfoObject class. + """ + def __init__(self, **attributes): + self.info = attributes + + def __getattr__(self, name): + try: + return self.info[name] + except KeyError: + raise AttributeError + + def get(self, key, default=None): + """Return the requested info. + """ + return self.info.get(key, default) + + +class Dataset(InfoObject): + """Dataset class. + """ + def __init__(self, data, **attributes): + InfoObject.__init__(self, **attributes) + self.data = data + + def __str__(self): + return str(self.data) + "\n" + str(self.info) + + def __repr__(self): + return repr(self.data) + "\n" + repr(self.info) + + def copy(self, copy_data=True): + """Return a copy of the data. + """ + if copy_data: + data = self.data.copy() + else: + data = self.data + return Dataset(data, **self.info) + + +class ProductList(object): + """Class for reading and storing product lists. + """ + def __init__(self, fname): + self.fname = fname + + tree = etree.parse(fname) + self._xml = tree.getroot() + self.prodlist = None + self.attrib = {} + self.vars = {} + self.aliases = {} + self.groups = [] + self.parse() + + def insert_vars(self): + """Variable replacement + """ + for item in self.prodlist.getiterator(): + for key in self.vars: + if key in item.attrib and item.attrib[key] in self.vars[key]: + item.set(key, self.vars[key][item.attrib[key]]) + + def check_groups(self): + """Check area groups. + """ + # create the "rest" group + last_group = [] + for area in self.prodlist: + assigned = False + for group in self.groups: + if len(area.attrib.keys()) == 0: + assigned = True + continue + if area.attrib["id"] in group.data: + assigned = True + break + if not assigned: + last_group.append(area.attrib["id"]) + if last_group: + self.groups.append(Dataset(last_group, id="_rest")) + + # replace area ids with actual xml area items + groups = [] + for group in self.groups: + new_group = Dataset([], **group.info) + for area_id in group.data: + assigned = False + for area in self.prodlist: + if area.attrib["id"] == area_id: + new_group.data.append(area) + assigned = True + break + if not assigned: + LOGGER.warning("Couldn't find area %s in product list", + area_id) + groups.append(new_group) + self.groups = groups + + def parse(self): + """Parse product list XML file. + """ + for item in self._xml: + if item.tag == "product_list": + self.prodlist = item + elif item.tag == "common": + for citem in item: + self.attrib[citem.tag] = citem.text + elif item.tag == "groups": + for group in item: + self.groups.append(Dataset(group.text.split(","), + **group.attrib)) + elif item.tag == "variables": + if item.attrib: + res = [os.environ[env] != val + for env, val in item.attrib.items()] + if any(res): + continue + for var in item: + self.vars.setdefault( + var.tag, {})[var.attrib["id"]] = var.text + elif item.tag == "aliases": + for alias in item: + self.aliases.setdefault( + alias.tag, + {})[alias.attrib["src"]] = alias.attrib['dst'] + self.insert_vars() + self.check_groups() + + +def get_root(fname): + '''Read XML file and return the root tree. + ''' + + tree = etree.parse(fname) + root = tree.getroot() + + return root + + +def parse_xml(tree, also_empty=False): + '''Parse the given XML file to dictionary. + ''' + xml_dict = {} + + # these tags will always be lists, if they are present and non-empty + listify = ['area', 'product', 'valid_satellite', 'invalid_satellite', + 'pattern', 'file_tag', 'directory'] + children = list(tree) + + if len(children) == 0: + try: + xml_dict = tree.text.strip() + except AttributeError: + pass + + for child in children: + new_val = parse_xml(child, also_empty=also_empty) + if len(new_val) == 0: + if also_empty: + xml_dict[child.tag] = '' + continue + else: + continue + + print "child.tag: ", child.tag + + if child.tag in xml_dict: + if not isinstance(xml_dict[child.tag], list): + xml_dict[child.tag] = [xml_dict[child.tag]] + xml_dict[child.tag].append(new_val) + else: + if len(new_val) > 0: + if child.tag in listify: + xml_dict[child.tag] = [new_val] + else: + xml_dict[child.tag] = new_val + + print xml_dict + return xml_dict + + +def get_filepattern_config(fname=None): + '''Retrieves the filepattern configuration file for trollstalker, + and returns the parsed XML as a dictionary. Optional argument + *fname* can be used to specify the file. If *fname* is None, the + systemwide file is read. + ''' + + if fname is None: + fname = os.path.realpath(__file__).split(os.path.sep)[:-2] + fname.append('etc') + fname.append('filepattern_config.xml') + fname = os.path.sep.join(fname) + + return parse_xml(get_root(fname), also_empty=True) From 64db8ee9a2e18374ae8f3862e0d7142bcbaf74ff Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Tue, 31 May 2016 08:20:57 +0200 Subject: [PATCH 04/27] Reading parts of config write sorted tiff --- mpop/satout/mitiff.py | 140 +++++++++++++----------------------------- 1 file changed, 41 insertions(+), 99 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index c3c6c814..684f7b1f 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -3,6 +3,7 @@ # Copyright (c) 2016. #from tifffile.tifffile import TIFF_TAGS from TiffImagePlugin import IMAGEDESCRIPTION +from libxml2mod import last # Author(s): @@ -37,89 +38,6 @@ logger = logging.getLogger(__name__) -#------------------------------------------------------------------------------- -# -# Read MiTIFF products config file. -# -#------------------------------------------------------------------------------- -def get_product_config(product_name, force_read=False): - """Read MiTIFF configuration entry for a given product name. - - :Parameters: - product_name : str - Name of MiTIFF product. - - :Arguments: - force_read : Boolean - Force re-reading config file. - - **Notes**: - * It will look for a *mitiff_products.cfg* in MPOP's - configuration directory defined by *PPP_CONFIG_DIR*. - * As an example, see *mitiff_products.cfg.template* in - MPOP's *etc* directory. - """ - return ProductConfigs()(product_name, force_read) - - -class _Singleton(type): - def __init__(cls, name_, bases_, dict_): - super(_Singleton, cls).__init__(name_, bases_, dict_) - cls.instance = None - - def __call__(cls, *args, **kwargs): - if cls.instance is None: - cls.instance = super(_Singleton, cls).__call__(*args, **kwargs) - return cls.instance - -class ProductConfigs(object): - __metaclass__ = _Singleton - - def __init__(self): - self.read_config() - - def __call__(self, product_name, force_read=False): - if force_read: - self.read_config() - return self._products[product_name] - - @property - def product_names(self): - return sorted(self._products.keys()) - - def read_config(self): - from ConfigParser import ConfigParser - - def _eval(val): - try: - return eval(val) - except: - return str(val) - - filename = self._find_a_config_file() - logger.info("Reading MiTIFF config file: '%s'" % filename) - - cfg = ConfigParser() - cfg.read(filename) - products = {} - for sec in cfg.sections(): - prd = {} - for key, val in cfg.items(sec): - prd[key] = _eval(val) - logger.debug("%s %s" % (key,prd[key])) - products[sec] = prd - self._products = products - - @staticmethod - def _find_a_config_file(): - name_ = 'mitiff_products.cfg' - home_ = os.path.dirname(os.path.abspath(__file__)) - penv_ = os.environ.get('PPP_CONFIG_DIR', '') - for fname_ in [os.path.join(x, name_) for x in (home_, penv_)]: - if os.path.isfile(fname_): - return fname_ - raise ValueError("Could not find a MiTIFF config file") - def save(scene, filename, product_name=None): """Saves the scene as a MITIFF file. This format can be read by DIANA @@ -132,21 +50,41 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N """ Write data to MITIFF file. """ print "Inside mitiff_writer ... " + import xml.etree.ElementTree as etree name_ = 'mitiff_products.xml' penv_ = os.environ.get('PPP_CONFIG_DIR', '') fname_ = os.path.join(penv_,name_) if os.path.isfile(fname_): - config = xml_read.parse_xml(xml_read.get_root(fname_)) + tree = etree.parse(fname_) + root = tree.getroot() + + #config = xml_read.parse_xml(xml_read.get_root(fname_)) #raise ValueError("Could not find a MiTIFF config file ", fname_) + + print "instrument_name: ", root_object.instrument_name + + product_list = root.findall("product[@instrument_name='%s']" % (root_object.instrument_name)) + #product_iter = root.iterfind("product[@instrument_name='%s']" % (root_object.instrument_name)) + product_list_length = len(product_list) - #if mitiff_product_name: - # options = deepcopy(get_product_config(mitiff_product_name)) - #else: - # print "No MiTIFF product name given. Skip Config. Sure you want this?" - # options = {} + if product_list_length == 0: + raise ValueError("No configuration mitiff product for this instrument is found: ", root_object.instrument_name) + elif product_list_length != 1: + raise ValueError("More than one configuration mitiff product for this instrument, '%s', found. Please check your config." % (root_object.instrument_name)) + + + for elem in product_list: + print "product_list : ", elem.tag, elem.attrib['instrument_name'] + product_channels = root.findall("product[@instrument_name='%s']/channel" % (root_object.instrument_name)) + for ch in product_channels: + print ch.items() + + #for child in root: + # print "loop: ", child + # print "loop2: ", child.tag, child.attrib #print root_object tifargs = dict() @@ -157,17 +95,21 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif = TIFF.open(filename, mode ='w') #print dir(tif) - #np.set_printoptions(threshold=np.nan) - - #tif.TIFFSetField(tif,IMAGEDESCRIPTION,"test") tif.SetField(IMAGEDESCRIPTION, "TEST") - for ch in root_object.channels: - print ch - print ch.info.get('units', 'None') - #print ch.data - data=ch.data.astype(np.uint8) - tif.write_image(data,) - + for ch in product_channels: + found_channel = False + print ch.attrib['name'] + for channels in root_object.channels: + if ch.attrib['name'] == channels.name: + data=channels.data.astype(np.uint8) + tif.write_image(data,) + found_channel = True + break + if not found_channel: + logger.debug("Could not find configured channel in read data set. Fill with empty.") + fill_channel = np.zeros(root_object.channels[0].data.shape,dtype=np.uint8) + tif.write_image(fill_channel) + tif.close if __name__ == '__main__': From 908b6fe31e4b904484a172a7579e22fb909d159d Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Tue, 31 May 2016 12:48:48 +0200 Subject: [PATCH 05/27] Writing scaled data to tiff --- mpop/satout/mitiff.py | 41 +++++++++++++++++++++++++++++++---------- 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 684f7b1f..853fd8ad 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -3,7 +3,8 @@ # Copyright (c) 2016. #from tifffile.tifffile import TIFF_TAGS from TiffImagePlugin import IMAGEDESCRIPTION -from libxml2mod import last +#from libxml2mod import last +#from aifc import data # Author(s): @@ -29,11 +30,13 @@ __revision__ = 0.1 +KELVIN_TO_CELSIUS = -273.15 + import os import numpy as np import logging from copy import deepcopy -import xml_read +#import xml_read #import xml.etree.ElementTree as etree logger = logging.getLogger(__name__) @@ -62,7 +65,7 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N #config = xml_read.parse_xml(xml_read.get_root(fname_)) #raise ValueError("Could not find a MiTIFF config file ", fname_) - print "instrument_name: ", root_object.instrument_name + #print "instrument_name: ", root_object.instrument_name product_list = root.findall("product[@instrument_name='%s']" % (root_object.instrument_name)) #product_iter = root.iterfind("product[@instrument_name='%s']" % (root_object.instrument_name)) @@ -82,10 +85,6 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N for ch in product_channels: print ch.items() - #for child in root: - # print "loop: ", child - # print "loop2: ", child.tag, child.attrib - #print root_object tifargs = dict() args={} @@ -95,21 +94,43 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif = TIFF.open(filename, mode ='w') #print dir(tif) + np.set_printoptions(threshold=np.nan) + tif.SetField(IMAGEDESCRIPTION, "TEST") for ch in product_channels: found_channel = False print ch.attrib['name'] for channels in root_object.channels: if ch.attrib['name'] == channels.name: - data=channels.data.astype(np.uint8) - tif.write_image(data,) + #Need to scale the data set to mitiff 0-255. 0 is no/missing/bad data. + logger.debug("min %f max %f value" % (float(ch.attrib['min-val']),float(ch.attrib['max-val']))) + reverse_offset = 0. + reverse_scale = 1. + if ch.attrib['calibration'] == "BT": + reverse_offset = 255. + reverse_scale = -1. + channels.data += KELVIN_TO_CELSIUS + #print channels.data + + logger.debug("Reverse offset: %f reverse scale: %f" % ( reverse_offset,reverse_scale)) + + _mask = channels.data.mask + _data = np.clip(channels.data, float(ch.attrib['min-val']),float(ch.attrib['max-val'])) + + data=reverse_offset + reverse_scale*(((_data-float(ch.attrib['min-val']))/(float(ch.attrib['max-val']) - float(ch.attrib['min-val'])))*255.) + + data[_mask] = 0 + + tif.write_image(data.astype(np.uint8),) found_channel = True break + if not found_channel: logger.debug("Could not find configured channel in read data set. Fill with empty.") fill_channel = np.zeros(root_object.channels[0].data.shape,dtype=np.uint8) tif.write_image(fill_channel) - + + tif.close if __name__ == '__main__': From 6b36ef20a68cb2f45ac3ce8a3a79efa8d705d3cc Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 1 Jun 2016 09:08:07 +0200 Subject: [PATCH 06/27] Working version with calibration --- mpop/satout/mitiff.py | 155 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 153 insertions(+), 2 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 853fd8ad..6312fbea 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -3,6 +3,7 @@ # Copyright (c) 2016. #from tifffile.tifffile import TIFF_TAGS from TiffImagePlugin import IMAGEDESCRIPTION +from pyparsing import Upcase #from libxml2mod import last #from aifc import data @@ -94,9 +95,15 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif = TIFF.open(filename, mode ='w') #print dir(tif) - np.set_printoptions(threshold=np.nan) + #np.set_printoptions(threshold=np.nan) - tif.SetField(IMAGEDESCRIPTION, "TEST") + #Need to generate information to go into the image description + #This is parsed by DIANA to be used when displaying the data + image_description = _make_image_description(root_object, product_channels) + + print str(image_description) + #print dtype(image_description) + tif.SetField(IMAGEDESCRIPTION, str(image_description)) for ch in product_channels: found_channel = False print ch.attrib['name'] @@ -133,6 +140,150 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif.close +def _make_image_description(ro, pc): + #generate image desdcription for mitiff. + """ + Satellite: NOAA 18 + Date and Time: 06:58 31/05-2016 + SatDir: 0 + Channels: 6 In this file: 1-VIS0.63 2-VIS0.86 3(3B)-IR3.7 4-IR10.8 5-IR11.5 6(3A)-VIS1.6 + Xsize: 4720 + Ysize: 5544 + Map projection: Stereographic + Proj string: +proj=stere +lon_0=0 +lat_0=90 +lat_ts=60 +ellps=WGS84 +towgs84=0,0,0 +units=km +x_0=2526000.000000 +y_0=5806000.000000 + TrueLat: 60 N + GridRot: 0 + Xunit:1000 m Yunit: 1000 m + NPX: 0.000000 NPY: 0.000000 + Ax: 1.000000 Ay: 1.000000 Bx: -2526.000000 By: -262.000000 + + Satellite: + Date and Time: + SatDir: 0 + Channels: In this file: + Xsize: + Ysize: + Map projection: Stereographic + Proj string: + TrueLat: 60 N + GridRot: 0 + Xunit:1000 m Yunit: 1000 m + NPX: 0.000000 NPY: 0.000000 + Ax: Ay: Bx: By: + + + if palette image write special palette + if normal channel write table calibration: + Table_calibration: , , [], , []\n\n + """ + + _image_description = '' + + _image_description += ' Satellite: ' + if ( ro.satname != None ): + _image_description += ro.satname.upper() + if ( ro.number != None ): + _image_description += ' '+str(ro.number) + + _image_description += '\n' + + _image_description += ' Date and Time: ' + _image_description += ro.time_slot.strftime("%H:%M %d/%m-%Y\n") + + _image_description += ' SatDir: 0\n' + + _image_description += ' Channels: ' + _image_description += str(len(pc)) + _image_description += ' In this file: ' + for ch in pc: + _image_description += ch.attrib['name-alias'] + if ch == pc[-1]: + _image_description += '\n' + else: + _image_description += ' ' + + _image_description += ' Xsize: ' + _image_description += str(ro.area.x_size) + '\n' + + _image_description += ' Ysize: ' + _image_description += str(ro.area.y_size) + '\n' + + _image_description += ' Map projection: Stereographic\n' + _image_description += ' Proj string: ' + ro.area.proj4_string + if not ("datum","towgs84") in ro.area.proj_dict: + _image_description += ' +towgs84=0,0,0' + + if not ("units") in ro.area.proj_dict: + _image_description += ' +units=km' + + + _image_description += ' +x_0=%.6f' % (-ro.area.area_extent[0]) + _image_description += ' +y_0=%.6f' % (-ro.area.area_extent[1]) + + _image_description += '\n' + _image_description += ' TrueLat: 60N\n' + _image_description += ' GridRot: 0\n' + + _image_description += ' Xunit:1000 m Yunit: 1000 m\n' + + _image_description += ' NPX: %.6f' % (0) + _image_description += ' NPY: %.6f' % (0) + '\n' + + _image_description += ' Ax: %.6f' % (ro.area.pixel_size_x/1000.) + _image_description += ' Ay: %.6f' % (ro.area.pixel_size_y/1000.) + #But this ads up to upper left corner of upper left pixel. + _image_description += ' Bx: %.6f' % (ro.area.area_extent[0]/1000.) #LL_x + _image_description += ' By: %.6f' % (ro.area.area_extent[3]/1000.) #UR_y + _image_description += '\n' + + logger.debug("Area extent: ",ro.area.area_extent) + + for ch in pc: + palette=False + #Make calibration. + if palette: + logging.debug("Do palette calibration") + else: + _image_description += 'Table_calibration: ' + if 'name-alias' in ch.keys(): + _image_description += ch.attrib['name-alias'] + else: + _image_description += ch.attrib['name'] + + _reverse_offset = 0.; + _reverse_scale = 1.; + + if ch.attrib['calibration'] == 'RADIANCE': + _image_description += ', Radiance, ' + _image_description += '[W/m²/µm/sr]' + _decimals = 8 + elif ch.attrib['calibration'] == 'BT': + _image_description += ', BT, ' + _image_description += '[°C]' + + _reverse_offset = 255.; + _reverse_scale = -1.; + _decimals = 2 + elif ch.attrib['calibration'] == 'REFLECTANCE': + _image_description += ', Reflectance(Albedo), ' + _image_description += '[%]' + _decimals = 2 + + else: + logger.warning("Unknown calib type. Must be Radiance, Reflectance or BT.") + + _image_description += ', 8, [ ' + for val in range(0,256): + #Comma separated list of values + #calib.append(boost::str(boost::format("%.8f ") % (prod_chan_it->min_val + (val * (prod_chan_it->max_val - prod_chan_it->min_val)) / 255.))); + _image_description += '{0:.{1}f} '.format((float(ch.attrib['min-val']) + ( (_reverse_offset + _reverse_scale*val) * ( float(ch.attrib['max-val']) - float(ch.attrib['min-val'])))/255.),_decimals) + + _image_description += ']\n\n' + + + return _image_description + + if __name__ == '__main__': from mpop.satellites import PolarFactory from mpop.projector import get_area_def From e3ddc7ee18392f499850e4772a0b9695156366cf Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 1 Jun 2016 09:09:05 +0200 Subject: [PATCH 07/27] Remove not needed file --- mpop/satout/xml_read.py | 230 ---------------------------------------- 1 file changed, 230 deletions(-) delete mode 100644 mpop/satout/xml_read.py diff --git a/mpop/satout/xml_read.py b/mpop/satout/xml_read.py deleted file mode 100644 index c729b90f..00000000 --- a/mpop/satout/xml_read.py +++ /dev/null @@ -1,230 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -# Copyright (c) 2014, 2015 - -# Author(s): - -# Panu Lahtinen -# Martin Raspaud - -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. - -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -'''XML reader for Trollduction system and product configuration files. -''' - -import xml.etree.ElementTree as etree -import os -import logging - -LOGGER = logging.getLogger(__name__) - - -class InfoObject(object): - """InfoObject class. - """ - def __init__(self, **attributes): - self.info = attributes - - def __getattr__(self, name): - try: - return self.info[name] - except KeyError: - raise AttributeError - - def get(self, key, default=None): - """Return the requested info. - """ - return self.info.get(key, default) - - -class Dataset(InfoObject): - """Dataset class. - """ - def __init__(self, data, **attributes): - InfoObject.__init__(self, **attributes) - self.data = data - - def __str__(self): - return str(self.data) + "\n" + str(self.info) - - def __repr__(self): - return repr(self.data) + "\n" + repr(self.info) - - def copy(self, copy_data=True): - """Return a copy of the data. - """ - if copy_data: - data = self.data.copy() - else: - data = self.data - return Dataset(data, **self.info) - - -class ProductList(object): - """Class for reading and storing product lists. - """ - def __init__(self, fname): - self.fname = fname - - tree = etree.parse(fname) - self._xml = tree.getroot() - self.prodlist = None - self.attrib = {} - self.vars = {} - self.aliases = {} - self.groups = [] - self.parse() - - def insert_vars(self): - """Variable replacement - """ - for item in self.prodlist.getiterator(): - for key in self.vars: - if key in item.attrib and item.attrib[key] in self.vars[key]: - item.set(key, self.vars[key][item.attrib[key]]) - - def check_groups(self): - """Check area groups. - """ - # create the "rest" group - last_group = [] - for area in self.prodlist: - assigned = False - for group in self.groups: - if len(area.attrib.keys()) == 0: - assigned = True - continue - if area.attrib["id"] in group.data: - assigned = True - break - if not assigned: - last_group.append(area.attrib["id"]) - if last_group: - self.groups.append(Dataset(last_group, id="_rest")) - - # replace area ids with actual xml area items - groups = [] - for group in self.groups: - new_group = Dataset([], **group.info) - for area_id in group.data: - assigned = False - for area in self.prodlist: - if area.attrib["id"] == area_id: - new_group.data.append(area) - assigned = True - break - if not assigned: - LOGGER.warning("Couldn't find area %s in product list", - area_id) - groups.append(new_group) - self.groups = groups - - def parse(self): - """Parse product list XML file. - """ - for item in self._xml: - if item.tag == "product_list": - self.prodlist = item - elif item.tag == "common": - for citem in item: - self.attrib[citem.tag] = citem.text - elif item.tag == "groups": - for group in item: - self.groups.append(Dataset(group.text.split(","), - **group.attrib)) - elif item.tag == "variables": - if item.attrib: - res = [os.environ[env] != val - for env, val in item.attrib.items()] - if any(res): - continue - for var in item: - self.vars.setdefault( - var.tag, {})[var.attrib["id"]] = var.text - elif item.tag == "aliases": - for alias in item: - self.aliases.setdefault( - alias.tag, - {})[alias.attrib["src"]] = alias.attrib['dst'] - self.insert_vars() - self.check_groups() - - -def get_root(fname): - '''Read XML file and return the root tree. - ''' - - tree = etree.parse(fname) - root = tree.getroot() - - return root - - -def parse_xml(tree, also_empty=False): - '''Parse the given XML file to dictionary. - ''' - xml_dict = {} - - # these tags will always be lists, if they are present and non-empty - listify = ['area', 'product', 'valid_satellite', 'invalid_satellite', - 'pattern', 'file_tag', 'directory'] - children = list(tree) - - if len(children) == 0: - try: - xml_dict = tree.text.strip() - except AttributeError: - pass - - for child in children: - new_val = parse_xml(child, also_empty=also_empty) - if len(new_val) == 0: - if also_empty: - xml_dict[child.tag] = '' - continue - else: - continue - - print "child.tag: ", child.tag - - if child.tag in xml_dict: - if not isinstance(xml_dict[child.tag], list): - xml_dict[child.tag] = [xml_dict[child.tag]] - xml_dict[child.tag].append(new_val) - else: - if len(new_val) > 0: - if child.tag in listify: - xml_dict[child.tag] = [new_val] - else: - xml_dict[child.tag] = new_val - - print xml_dict - return xml_dict - - -def get_filepattern_config(fname=None): - '''Retrieves the filepattern configuration file for trollstalker, - and returns the parsed XML as a dictionary. Optional argument - *fname* can be used to specify the file. If *fname* is None, the - systemwide file is read. - ''' - - if fname is None: - fname = os.path.realpath(__file__).split(os.path.sep)[:-2] - fname.append('etc') - fname.append('filepattern_config.xml') - fname = os.path.sep.join(fname) - - return parse_xml(get_root(fname), also_empty=True) From 9d26b8065422307bd9a02a68573db1f9372878a6 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 1 Jun 2016 10:19:30 +0200 Subject: [PATCH 08/27] Bumped version Added IOError on not found xml file --- mpop/satout/mitiff.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 6312fbea..2bbb3a54 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -29,7 +29,7 @@ """mpop mitiff writer interface. """ -__revision__ = 0.1 +__revision__ = 0.2 KELVIN_TO_CELSIUS = -273.15 @@ -42,17 +42,17 @@ logger = logging.getLogger(__name__) -def save(scene, filename, product_name=None): +def save(scene, filename): """Saves the scene as a MITIFF file. This format can be read by DIANA """ - print "inside save ... " - return mitiff_writer(filename, scene, mitiff_product_name=product_name) + logger.debug("inside save ... ") + return mitiff_writer(filename, scene) -def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=None): +def mitiff_writer(filename, root_object, compression=True): """ Write data to MITIFF file. """ - print "Inside mitiff_writer ... " + logger.debug("Inside mitiff_writer ... ") import xml.etree.ElementTree as etree @@ -62,10 +62,9 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N if os.path.isfile(fname_): tree = etree.parse(fname_) root = tree.getroot() - - #config = xml_read.parse_xml(xml_read.get_root(fname_)) - #raise ValueError("Could not find a MiTIFF config file ", fname_) - + else: + raise IOError("Could not find %s. Check if this file is available at this location." % (fname_)) + #print "instrument_name: ", root_object.instrument_name product_list = root.findall("product[@instrument_name='%s']" % (root_object.instrument_name)) @@ -272,6 +271,9 @@ def _make_image_description(ro, pc): else: logger.warning("Unknown calib type. Must be Radiance, Reflectance or BT.") + #How to format string by passing the format + #http://stackoverflow.com/questions/1598579/rounding-decimals-with-new-python-format-function + _image_description += ', 8, [ ' for val in range(0,256): #Comma separated list of values @@ -306,7 +308,3 @@ def _make_image_description(ro, pc): save(local_data,"test.mitiff",product_name="AVHRR") print "Complete." - - - - \ No newline at end of file From 53e900a0802a2e538053988b00084efa6c2d82c9 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:40:53 +0000 Subject: [PATCH 09/27] mpop/satout/mitiff.py some more debug info --- mpop/satout/mitiff.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 2bbb3a54..197f4bf0 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -42,7 +42,7 @@ logger = logging.getLogger(__name__) -def save(scene, filename): +def save(scene, filename, **options): """Saves the scene as a MITIFF file. This format can be read by DIANA """ @@ -235,13 +235,13 @@ def _make_image_description(ro, pc): _image_description += ' By: %.6f' % (ro.area.area_extent[3]/1000.) #UR_y _image_description += '\n' - logger.debug("Area extent: ",ro.area.area_extent) + logger.debug("Area extent: {}".format(ro.area.area_extent)) for ch in pc: palette=False #Make calibration. if palette: - logging.debug("Do palette calibration") + logging.debug("Do palette calibration. Not implemented yet.") else: _image_description += 'Table_calibration: ' if 'name-alias' in ch.keys(): From f2eeb9113ee36a6bf4a60e371251a07e4d22ce25 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:41:33 +0000 Subject: [PATCH 10/27] mpop/imageo/geo_image.py more info --- mpop/imageo/geo_image.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/mpop/imageo/geo_image.py b/mpop/imageo/geo_image.py index c1157b03..7b1316da 100644 --- a/mpop/imageo/geo_image.py +++ b/mpop/imageo/geo_image.py @@ -96,6 +96,9 @@ def save(self, filename, compression=6, if fformat.lower() in ('tif', 'tiff'): return self.geotiff_save(filename, compression, tags, gdal_options, blocksize, **kwargs) + elif fformat.lower() in ('mitiff'): + logger.info("Writing to mitiff ....") + try: # Let image.pil_save it ? Image.save(self, filename, compression, fformat=fformat) From d9ae21336804243f1f3120f386a17d82c06c6e54 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:44:29 +0000 Subject: [PATCH 11/27] mpop/scene.py changed save definition. This could possible break the code. This needs to be confirmed --- mpop/scene.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mpop/scene.py b/mpop/scene.py index 9d993ef1..3fa046fc 100644 --- a/mpop/scene.py +++ b/mpop/scene.py @@ -499,14 +499,16 @@ def load(self, channels=None, load_again=False, area_extent=None, **kwargs): self.channels_to_load = set() - def save(self, filename, to_format="netcdf4", **options): + def save(self, filename, fformat="netcdf4", compression=None, **options): """Saves the current scene into a file of format *to_format*. Supported formats are: - *netcdf4*: NetCDF4 with CF conventions. """ - + to_format = fformat writer = "satout." + to_format + LOG.debug("writer module: %s" %(writer)) + try: writer_module = __import__(writer, globals(), locals(), ["save"]) except ImportError, err: From 405f153fd2f909549fc3dfdfb6500bf1eedd7825 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:46:47 +0000 Subject: [PATCH 12/27] mpop/scene.py added info in mitiff save format --- mpop/scene.py | 1 + 1 file changed, 1 insertion(+) diff --git a/mpop/scene.py b/mpop/scene.py index 3fa046fc..c14127f8 100644 --- a/mpop/scene.py +++ b/mpop/scene.py @@ -504,6 +504,7 @@ def save(self, filename, fformat="netcdf4", compression=None, **options): formats are: - *netcdf4*: NetCDF4 with CF conventions. + - *mitiff*: MITIFF format, TIFF file with image description tag holding the geolocation info. Used by DIANA at Norwegian MET """ to_format = fformat writer = "satout." + to_format From a3be102231d85ed79b6bc2cf02929f6108c8f728 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 25 May 2016 07:35:47 +0000 Subject: [PATCH 13/27] mpop/satout/mitiff.py first draft --- mpop/satout/mitiff.py | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 mpop/satout/mitiff.py diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py new file mode 100644 index 00000000..0d9cf108 --- /dev/null +++ b/mpop/satout/mitiff.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2016. + +# Author(s): + +# Trygve Aspenes + +# This file is part of mpop. + +# mpop is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# mpop is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License along with +# mpop. If not, see . + + +"""mpop mitiff writer interface. +""" + +__revision__ = 0.1 + +import numpy as np +import logging + +logger = logging.getLogger(__name__) + +def save(scene, filename, compression=True, dtype=np.int16, band_axis=2): + """Saves the scene as a MITIFF file. This is local for METNO, but can be read by DIANA + + """ + return mitiff_writer(filename, CFScene(scene, dtype, band_axis), compression=compression) + + +def mitiff__writer(filename, root_object, compression=True): + """ Write data to MITIFF file. """ + From 5b2fe143106896c3a8d54ab53ab7475e1e58a8cd Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Mon, 30 May 2016 09:52:21 +0200 Subject: [PATCH 14/27] Write files ok. Added reading of xml config file, work in progress remove reading of cfg style config --- mpop/satout/mitiff.py | 165 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 160 insertions(+), 5 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 0d9cf108..c3c6c814 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -1,6 +1,8 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # Copyright (c) 2016. +#from tifffile.tifffile import TIFF_TAGS +from TiffImagePlugin import IMAGEDESCRIPTION # Author(s): @@ -26,18 +28,171 @@ __revision__ = 0.1 +import os import numpy as np import logging +from copy import deepcopy +import xml_read +#import xml.etree.ElementTree as etree logger = logging.getLogger(__name__) -def save(scene, filename, compression=True, dtype=np.int16, band_axis=2): - """Saves the scene as a MITIFF file. This is local for METNO, but can be read by DIANA +#------------------------------------------------------------------------------- +# +# Read MiTIFF products config file. +# +#------------------------------------------------------------------------------- +def get_product_config(product_name, force_read=False): + """Read MiTIFF configuration entry for a given product name. + + :Parameters: + product_name : str + Name of MiTIFF product. + + :Arguments: + force_read : Boolean + Force re-reading config file. + + **Notes**: + * It will look for a *mitiff_products.cfg* in MPOP's + configuration directory defined by *PPP_CONFIG_DIR*. + * As an example, see *mitiff_products.cfg.template* in + MPOP's *etc* directory. + """ + return ProductConfigs()(product_name, force_read) + + +class _Singleton(type): + def __init__(cls, name_, bases_, dict_): + super(_Singleton, cls).__init__(name_, bases_, dict_) + cls.instance = None + + def __call__(cls, *args, **kwargs): + if cls.instance is None: + cls.instance = super(_Singleton, cls).__call__(*args, **kwargs) + return cls.instance + +class ProductConfigs(object): + __metaclass__ = _Singleton + + def __init__(self): + self.read_config() + + def __call__(self, product_name, force_read=False): + if force_read: + self.read_config() + return self._products[product_name] + + @property + def product_names(self): + return sorted(self._products.keys()) + + def read_config(self): + from ConfigParser import ConfigParser + + def _eval(val): + try: + return eval(val) + except: + return str(val) + + filename = self._find_a_config_file() + logger.info("Reading MiTIFF config file: '%s'" % filename) + + cfg = ConfigParser() + cfg.read(filename) + products = {} + for sec in cfg.sections(): + prd = {} + for key, val in cfg.items(sec): + prd[key] = _eval(val) + logger.debug("%s %s" % (key,prd[key])) + products[sec] = prd + self._products = products + + @staticmethod + def _find_a_config_file(): + name_ = 'mitiff_products.cfg' + home_ = os.path.dirname(os.path.abspath(__file__)) + penv_ = os.environ.get('PPP_CONFIG_DIR', '') + for fname_ in [os.path.join(x, name_) for x in (home_, penv_)]: + if os.path.isfile(fname_): + return fname_ + raise ValueError("Could not find a MiTIFF config file") + +def save(scene, filename, product_name=None): + """Saves the scene as a MITIFF file. This format can be read by DIANA """ - return mitiff_writer(filename, CFScene(scene, dtype, band_axis), compression=compression) + print "inside save ... " + return mitiff_writer(filename, scene, mitiff_product_name=product_name) -def mitiff__writer(filename, root_object, compression=True): +def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=None): """ Write data to MITIFF file. """ - + print "Inside mitiff_writer ... " + + + name_ = 'mitiff_products.xml' + penv_ = os.environ.get('PPP_CONFIG_DIR', '') + fname_ = os.path.join(penv_,name_) + if os.path.isfile(fname_): + config = xml_read.parse_xml(xml_read.get_root(fname_)) + #raise ValueError("Could not find a MiTIFF config file ", fname_) + + + #if mitiff_product_name: + # options = deepcopy(get_product_config(mitiff_product_name)) + #else: + # print "No MiTIFF product name given. Skip Config. Sure you want this?" + # options = {} + + + #print root_object + tifargs = dict() + args={} + + from libtiff import TIFF + + tif = TIFF.open(filename, mode ='w') + #print dir(tif) + + #np.set_printoptions(threshold=np.nan) + + #tif.TIFFSetField(tif,IMAGEDESCRIPTION,"test") + tif.SetField(IMAGEDESCRIPTION, "TEST") + for ch in root_object.channels: + print ch + print ch.info.get('units', 'None') + #print ch.data + data=ch.data.astype(np.uint8) + tif.write_image(data,) + + tif.close + +if __name__ == '__main__': + from mpop.satellites import PolarFactory + from mpop.projector import get_area_def + import datetime + + logging.basicConfig(level=logging.DEBUG) + + print "Test saving mitiff." + orbit = "37582" + time_slot = datetime.datetime(2016, 5, 24, 15, 25) + global_data = PolarFactory.create_scene("noaa","19","avhrr/3",time_slot, orbit) + global_data.load() + + from pprint import pprint + #pprint(vars(global_data)) + + local_data = global_data.project("eurol_metno", mode="nearest") + + pprint(vars(local_data)) + + save(local_data,"test.mitiff",product_name="AVHRR") + print "Complete." + + + + \ No newline at end of file From f373addb97c44cc0c606789bc8049cea6882277f Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Mon, 30 May 2016 10:05:03 +0200 Subject: [PATCH 15/27] xml_read taken from trollduction and must be adapted --- mpop/satout/xml_read.py | 230 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 230 insertions(+) create mode 100644 mpop/satout/xml_read.py diff --git a/mpop/satout/xml_read.py b/mpop/satout/xml_read.py new file mode 100644 index 00000000..c729b90f --- /dev/null +++ b/mpop/satout/xml_read.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2014, 2015 + +# Author(s): + +# Panu Lahtinen +# Martin Raspaud + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +'''XML reader for Trollduction system and product configuration files. +''' + +import xml.etree.ElementTree as etree +import os +import logging + +LOGGER = logging.getLogger(__name__) + + +class InfoObject(object): + """InfoObject class. + """ + def __init__(self, **attributes): + self.info = attributes + + def __getattr__(self, name): + try: + return self.info[name] + except KeyError: + raise AttributeError + + def get(self, key, default=None): + """Return the requested info. + """ + return self.info.get(key, default) + + +class Dataset(InfoObject): + """Dataset class. + """ + def __init__(self, data, **attributes): + InfoObject.__init__(self, **attributes) + self.data = data + + def __str__(self): + return str(self.data) + "\n" + str(self.info) + + def __repr__(self): + return repr(self.data) + "\n" + repr(self.info) + + def copy(self, copy_data=True): + """Return a copy of the data. + """ + if copy_data: + data = self.data.copy() + else: + data = self.data + return Dataset(data, **self.info) + + +class ProductList(object): + """Class for reading and storing product lists. + """ + def __init__(self, fname): + self.fname = fname + + tree = etree.parse(fname) + self._xml = tree.getroot() + self.prodlist = None + self.attrib = {} + self.vars = {} + self.aliases = {} + self.groups = [] + self.parse() + + def insert_vars(self): + """Variable replacement + """ + for item in self.prodlist.getiterator(): + for key in self.vars: + if key in item.attrib and item.attrib[key] in self.vars[key]: + item.set(key, self.vars[key][item.attrib[key]]) + + def check_groups(self): + """Check area groups. + """ + # create the "rest" group + last_group = [] + for area in self.prodlist: + assigned = False + for group in self.groups: + if len(area.attrib.keys()) == 0: + assigned = True + continue + if area.attrib["id"] in group.data: + assigned = True + break + if not assigned: + last_group.append(area.attrib["id"]) + if last_group: + self.groups.append(Dataset(last_group, id="_rest")) + + # replace area ids with actual xml area items + groups = [] + for group in self.groups: + new_group = Dataset([], **group.info) + for area_id in group.data: + assigned = False + for area in self.prodlist: + if area.attrib["id"] == area_id: + new_group.data.append(area) + assigned = True + break + if not assigned: + LOGGER.warning("Couldn't find area %s in product list", + area_id) + groups.append(new_group) + self.groups = groups + + def parse(self): + """Parse product list XML file. + """ + for item in self._xml: + if item.tag == "product_list": + self.prodlist = item + elif item.tag == "common": + for citem in item: + self.attrib[citem.tag] = citem.text + elif item.tag == "groups": + for group in item: + self.groups.append(Dataset(group.text.split(","), + **group.attrib)) + elif item.tag == "variables": + if item.attrib: + res = [os.environ[env] != val + for env, val in item.attrib.items()] + if any(res): + continue + for var in item: + self.vars.setdefault( + var.tag, {})[var.attrib["id"]] = var.text + elif item.tag == "aliases": + for alias in item: + self.aliases.setdefault( + alias.tag, + {})[alias.attrib["src"]] = alias.attrib['dst'] + self.insert_vars() + self.check_groups() + + +def get_root(fname): + '''Read XML file and return the root tree. + ''' + + tree = etree.parse(fname) + root = tree.getroot() + + return root + + +def parse_xml(tree, also_empty=False): + '''Parse the given XML file to dictionary. + ''' + xml_dict = {} + + # these tags will always be lists, if they are present and non-empty + listify = ['area', 'product', 'valid_satellite', 'invalid_satellite', + 'pattern', 'file_tag', 'directory'] + children = list(tree) + + if len(children) == 0: + try: + xml_dict = tree.text.strip() + except AttributeError: + pass + + for child in children: + new_val = parse_xml(child, also_empty=also_empty) + if len(new_val) == 0: + if also_empty: + xml_dict[child.tag] = '' + continue + else: + continue + + print "child.tag: ", child.tag + + if child.tag in xml_dict: + if not isinstance(xml_dict[child.tag], list): + xml_dict[child.tag] = [xml_dict[child.tag]] + xml_dict[child.tag].append(new_val) + else: + if len(new_val) > 0: + if child.tag in listify: + xml_dict[child.tag] = [new_val] + else: + xml_dict[child.tag] = new_val + + print xml_dict + return xml_dict + + +def get_filepattern_config(fname=None): + '''Retrieves the filepattern configuration file for trollstalker, + and returns the parsed XML as a dictionary. Optional argument + *fname* can be used to specify the file. If *fname* is None, the + systemwide file is read. + ''' + + if fname is None: + fname = os.path.realpath(__file__).split(os.path.sep)[:-2] + fname.append('etc') + fname.append('filepattern_config.xml') + fname = os.path.sep.join(fname) + + return parse_xml(get_root(fname), also_empty=True) From 430324a5c95939f6d53da9ead080ba8dd6de7f27 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Tue, 31 May 2016 08:20:57 +0200 Subject: [PATCH 16/27] Reading parts of config write sorted tiff --- mpop/satout/mitiff.py | 140 +++++++++++++----------------------------- 1 file changed, 41 insertions(+), 99 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index c3c6c814..684f7b1f 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -3,6 +3,7 @@ # Copyright (c) 2016. #from tifffile.tifffile import TIFF_TAGS from TiffImagePlugin import IMAGEDESCRIPTION +from libxml2mod import last # Author(s): @@ -37,89 +38,6 @@ logger = logging.getLogger(__name__) -#------------------------------------------------------------------------------- -# -# Read MiTIFF products config file. -# -#------------------------------------------------------------------------------- -def get_product_config(product_name, force_read=False): - """Read MiTIFF configuration entry for a given product name. - - :Parameters: - product_name : str - Name of MiTIFF product. - - :Arguments: - force_read : Boolean - Force re-reading config file. - - **Notes**: - * It will look for a *mitiff_products.cfg* in MPOP's - configuration directory defined by *PPP_CONFIG_DIR*. - * As an example, see *mitiff_products.cfg.template* in - MPOP's *etc* directory. - """ - return ProductConfigs()(product_name, force_read) - - -class _Singleton(type): - def __init__(cls, name_, bases_, dict_): - super(_Singleton, cls).__init__(name_, bases_, dict_) - cls.instance = None - - def __call__(cls, *args, **kwargs): - if cls.instance is None: - cls.instance = super(_Singleton, cls).__call__(*args, **kwargs) - return cls.instance - -class ProductConfigs(object): - __metaclass__ = _Singleton - - def __init__(self): - self.read_config() - - def __call__(self, product_name, force_read=False): - if force_read: - self.read_config() - return self._products[product_name] - - @property - def product_names(self): - return sorted(self._products.keys()) - - def read_config(self): - from ConfigParser import ConfigParser - - def _eval(val): - try: - return eval(val) - except: - return str(val) - - filename = self._find_a_config_file() - logger.info("Reading MiTIFF config file: '%s'" % filename) - - cfg = ConfigParser() - cfg.read(filename) - products = {} - for sec in cfg.sections(): - prd = {} - for key, val in cfg.items(sec): - prd[key] = _eval(val) - logger.debug("%s %s" % (key,prd[key])) - products[sec] = prd - self._products = products - - @staticmethod - def _find_a_config_file(): - name_ = 'mitiff_products.cfg' - home_ = os.path.dirname(os.path.abspath(__file__)) - penv_ = os.environ.get('PPP_CONFIG_DIR', '') - for fname_ in [os.path.join(x, name_) for x in (home_, penv_)]: - if os.path.isfile(fname_): - return fname_ - raise ValueError("Could not find a MiTIFF config file") - def save(scene, filename, product_name=None): """Saves the scene as a MITIFF file. This format can be read by DIANA @@ -132,21 +50,41 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N """ Write data to MITIFF file. """ print "Inside mitiff_writer ... " + import xml.etree.ElementTree as etree name_ = 'mitiff_products.xml' penv_ = os.environ.get('PPP_CONFIG_DIR', '') fname_ = os.path.join(penv_,name_) if os.path.isfile(fname_): - config = xml_read.parse_xml(xml_read.get_root(fname_)) + tree = etree.parse(fname_) + root = tree.getroot() + + #config = xml_read.parse_xml(xml_read.get_root(fname_)) #raise ValueError("Could not find a MiTIFF config file ", fname_) + + print "instrument_name: ", root_object.instrument_name + + product_list = root.findall("product[@instrument_name='%s']" % (root_object.instrument_name)) + #product_iter = root.iterfind("product[@instrument_name='%s']" % (root_object.instrument_name)) + product_list_length = len(product_list) - #if mitiff_product_name: - # options = deepcopy(get_product_config(mitiff_product_name)) - #else: - # print "No MiTIFF product name given. Skip Config. Sure you want this?" - # options = {} + if product_list_length == 0: + raise ValueError("No configuration mitiff product for this instrument is found: ", root_object.instrument_name) + elif product_list_length != 1: + raise ValueError("More than one configuration mitiff product for this instrument, '%s', found. Please check your config." % (root_object.instrument_name)) + + + for elem in product_list: + print "product_list : ", elem.tag, elem.attrib['instrument_name'] + product_channels = root.findall("product[@instrument_name='%s']/channel" % (root_object.instrument_name)) + for ch in product_channels: + print ch.items() + + #for child in root: + # print "loop: ", child + # print "loop2: ", child.tag, child.attrib #print root_object tifargs = dict() @@ -157,17 +95,21 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif = TIFF.open(filename, mode ='w') #print dir(tif) - #np.set_printoptions(threshold=np.nan) - - #tif.TIFFSetField(tif,IMAGEDESCRIPTION,"test") tif.SetField(IMAGEDESCRIPTION, "TEST") - for ch in root_object.channels: - print ch - print ch.info.get('units', 'None') - #print ch.data - data=ch.data.astype(np.uint8) - tif.write_image(data,) - + for ch in product_channels: + found_channel = False + print ch.attrib['name'] + for channels in root_object.channels: + if ch.attrib['name'] == channels.name: + data=channels.data.astype(np.uint8) + tif.write_image(data,) + found_channel = True + break + if not found_channel: + logger.debug("Could not find configured channel in read data set. Fill with empty.") + fill_channel = np.zeros(root_object.channels[0].data.shape,dtype=np.uint8) + tif.write_image(fill_channel) + tif.close if __name__ == '__main__': From 3b5be55c611e7b53933bcfa1f36b09d47cfdcec0 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Tue, 31 May 2016 12:48:48 +0200 Subject: [PATCH 17/27] Writing scaled data to tiff --- mpop/satout/mitiff.py | 41 +++++++++++++++++++++++++++++++---------- 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 684f7b1f..853fd8ad 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -3,7 +3,8 @@ # Copyright (c) 2016. #from tifffile.tifffile import TIFF_TAGS from TiffImagePlugin import IMAGEDESCRIPTION -from libxml2mod import last +#from libxml2mod import last +#from aifc import data # Author(s): @@ -29,11 +30,13 @@ __revision__ = 0.1 +KELVIN_TO_CELSIUS = -273.15 + import os import numpy as np import logging from copy import deepcopy -import xml_read +#import xml_read #import xml.etree.ElementTree as etree logger = logging.getLogger(__name__) @@ -62,7 +65,7 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N #config = xml_read.parse_xml(xml_read.get_root(fname_)) #raise ValueError("Could not find a MiTIFF config file ", fname_) - print "instrument_name: ", root_object.instrument_name + #print "instrument_name: ", root_object.instrument_name product_list = root.findall("product[@instrument_name='%s']" % (root_object.instrument_name)) #product_iter = root.iterfind("product[@instrument_name='%s']" % (root_object.instrument_name)) @@ -82,10 +85,6 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N for ch in product_channels: print ch.items() - #for child in root: - # print "loop: ", child - # print "loop2: ", child.tag, child.attrib - #print root_object tifargs = dict() args={} @@ -95,21 +94,43 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif = TIFF.open(filename, mode ='w') #print dir(tif) + np.set_printoptions(threshold=np.nan) + tif.SetField(IMAGEDESCRIPTION, "TEST") for ch in product_channels: found_channel = False print ch.attrib['name'] for channels in root_object.channels: if ch.attrib['name'] == channels.name: - data=channels.data.astype(np.uint8) - tif.write_image(data,) + #Need to scale the data set to mitiff 0-255. 0 is no/missing/bad data. + logger.debug("min %f max %f value" % (float(ch.attrib['min-val']),float(ch.attrib['max-val']))) + reverse_offset = 0. + reverse_scale = 1. + if ch.attrib['calibration'] == "BT": + reverse_offset = 255. + reverse_scale = -1. + channels.data += KELVIN_TO_CELSIUS + #print channels.data + + logger.debug("Reverse offset: %f reverse scale: %f" % ( reverse_offset,reverse_scale)) + + _mask = channels.data.mask + _data = np.clip(channels.data, float(ch.attrib['min-val']),float(ch.attrib['max-val'])) + + data=reverse_offset + reverse_scale*(((_data-float(ch.attrib['min-val']))/(float(ch.attrib['max-val']) - float(ch.attrib['min-val'])))*255.) + + data[_mask] = 0 + + tif.write_image(data.astype(np.uint8),) found_channel = True break + if not found_channel: logger.debug("Could not find configured channel in read data set. Fill with empty.") fill_channel = np.zeros(root_object.channels[0].data.shape,dtype=np.uint8) tif.write_image(fill_channel) - + + tif.close if __name__ == '__main__': From 98fe975df3592d6cdda19d7698ff10b981216e57 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 1 Jun 2016 09:08:07 +0200 Subject: [PATCH 18/27] Working version with calibration --- mpop/satout/mitiff.py | 155 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 153 insertions(+), 2 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 853fd8ad..6312fbea 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -3,6 +3,7 @@ # Copyright (c) 2016. #from tifffile.tifffile import TIFF_TAGS from TiffImagePlugin import IMAGEDESCRIPTION +from pyparsing import Upcase #from libxml2mod import last #from aifc import data @@ -94,9 +95,15 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif = TIFF.open(filename, mode ='w') #print dir(tif) - np.set_printoptions(threshold=np.nan) + #np.set_printoptions(threshold=np.nan) - tif.SetField(IMAGEDESCRIPTION, "TEST") + #Need to generate information to go into the image description + #This is parsed by DIANA to be used when displaying the data + image_description = _make_image_description(root_object, product_channels) + + print str(image_description) + #print dtype(image_description) + tif.SetField(IMAGEDESCRIPTION, str(image_description)) for ch in product_channels: found_channel = False print ch.attrib['name'] @@ -133,6 +140,150 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N tif.close +def _make_image_description(ro, pc): + #generate image desdcription for mitiff. + """ + Satellite: NOAA 18 + Date and Time: 06:58 31/05-2016 + SatDir: 0 + Channels: 6 In this file: 1-VIS0.63 2-VIS0.86 3(3B)-IR3.7 4-IR10.8 5-IR11.5 6(3A)-VIS1.6 + Xsize: 4720 + Ysize: 5544 + Map projection: Stereographic + Proj string: +proj=stere +lon_0=0 +lat_0=90 +lat_ts=60 +ellps=WGS84 +towgs84=0,0,0 +units=km +x_0=2526000.000000 +y_0=5806000.000000 + TrueLat: 60 N + GridRot: 0 + Xunit:1000 m Yunit: 1000 m + NPX: 0.000000 NPY: 0.000000 + Ax: 1.000000 Ay: 1.000000 Bx: -2526.000000 By: -262.000000 + + Satellite: + Date and Time: + SatDir: 0 + Channels: In this file: + Xsize: + Ysize: + Map projection: Stereographic + Proj string: + TrueLat: 60 N + GridRot: 0 + Xunit:1000 m Yunit: 1000 m + NPX: 0.000000 NPY: 0.000000 + Ax: Ay: Bx: By: + + + if palette image write special palette + if normal channel write table calibration: + Table_calibration: , , [], , []\n\n + """ + + _image_description = '' + + _image_description += ' Satellite: ' + if ( ro.satname != None ): + _image_description += ro.satname.upper() + if ( ro.number != None ): + _image_description += ' '+str(ro.number) + + _image_description += '\n' + + _image_description += ' Date and Time: ' + _image_description += ro.time_slot.strftime("%H:%M %d/%m-%Y\n") + + _image_description += ' SatDir: 0\n' + + _image_description += ' Channels: ' + _image_description += str(len(pc)) + _image_description += ' In this file: ' + for ch in pc: + _image_description += ch.attrib['name-alias'] + if ch == pc[-1]: + _image_description += '\n' + else: + _image_description += ' ' + + _image_description += ' Xsize: ' + _image_description += str(ro.area.x_size) + '\n' + + _image_description += ' Ysize: ' + _image_description += str(ro.area.y_size) + '\n' + + _image_description += ' Map projection: Stereographic\n' + _image_description += ' Proj string: ' + ro.area.proj4_string + if not ("datum","towgs84") in ro.area.proj_dict: + _image_description += ' +towgs84=0,0,0' + + if not ("units") in ro.area.proj_dict: + _image_description += ' +units=km' + + + _image_description += ' +x_0=%.6f' % (-ro.area.area_extent[0]) + _image_description += ' +y_0=%.6f' % (-ro.area.area_extent[1]) + + _image_description += '\n' + _image_description += ' TrueLat: 60N\n' + _image_description += ' GridRot: 0\n' + + _image_description += ' Xunit:1000 m Yunit: 1000 m\n' + + _image_description += ' NPX: %.6f' % (0) + _image_description += ' NPY: %.6f' % (0) + '\n' + + _image_description += ' Ax: %.6f' % (ro.area.pixel_size_x/1000.) + _image_description += ' Ay: %.6f' % (ro.area.pixel_size_y/1000.) + #But this ads up to upper left corner of upper left pixel. + _image_description += ' Bx: %.6f' % (ro.area.area_extent[0]/1000.) #LL_x + _image_description += ' By: %.6f' % (ro.area.area_extent[3]/1000.) #UR_y + _image_description += '\n' + + logger.debug("Area extent: ",ro.area.area_extent) + + for ch in pc: + palette=False + #Make calibration. + if palette: + logging.debug("Do palette calibration") + else: + _image_description += 'Table_calibration: ' + if 'name-alias' in ch.keys(): + _image_description += ch.attrib['name-alias'] + else: + _image_description += ch.attrib['name'] + + _reverse_offset = 0.; + _reverse_scale = 1.; + + if ch.attrib['calibration'] == 'RADIANCE': + _image_description += ', Radiance, ' + _image_description += '[W/m²/µm/sr]' + _decimals = 8 + elif ch.attrib['calibration'] == 'BT': + _image_description += ', BT, ' + _image_description += '[°C]' + + _reverse_offset = 255.; + _reverse_scale = -1.; + _decimals = 2 + elif ch.attrib['calibration'] == 'REFLECTANCE': + _image_description += ', Reflectance(Albedo), ' + _image_description += '[%]' + _decimals = 2 + + else: + logger.warning("Unknown calib type. Must be Radiance, Reflectance or BT.") + + _image_description += ', 8, [ ' + for val in range(0,256): + #Comma separated list of values + #calib.append(boost::str(boost::format("%.8f ") % (prod_chan_it->min_val + (val * (prod_chan_it->max_val - prod_chan_it->min_val)) / 255.))); + _image_description += '{0:.{1}f} '.format((float(ch.attrib['min-val']) + ( (_reverse_offset + _reverse_scale*val) * ( float(ch.attrib['max-val']) - float(ch.attrib['min-val'])))/255.),_decimals) + + _image_description += ']\n\n' + + + return _image_description + + if __name__ == '__main__': from mpop.satellites import PolarFactory from mpop.projector import get_area_def From 8d21cc765c9fe6bf5cc6bf04619680eaf4dd904e Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 1 Jun 2016 09:09:05 +0200 Subject: [PATCH 19/27] Remove not needed file --- mpop/satout/xml_read.py | 230 ---------------------------------------- 1 file changed, 230 deletions(-) delete mode 100644 mpop/satout/xml_read.py diff --git a/mpop/satout/xml_read.py b/mpop/satout/xml_read.py deleted file mode 100644 index c729b90f..00000000 --- a/mpop/satout/xml_read.py +++ /dev/null @@ -1,230 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -# Copyright (c) 2014, 2015 - -# Author(s): - -# Panu Lahtinen -# Martin Raspaud - -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. - -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. - -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -'''XML reader for Trollduction system and product configuration files. -''' - -import xml.etree.ElementTree as etree -import os -import logging - -LOGGER = logging.getLogger(__name__) - - -class InfoObject(object): - """InfoObject class. - """ - def __init__(self, **attributes): - self.info = attributes - - def __getattr__(self, name): - try: - return self.info[name] - except KeyError: - raise AttributeError - - def get(self, key, default=None): - """Return the requested info. - """ - return self.info.get(key, default) - - -class Dataset(InfoObject): - """Dataset class. - """ - def __init__(self, data, **attributes): - InfoObject.__init__(self, **attributes) - self.data = data - - def __str__(self): - return str(self.data) + "\n" + str(self.info) - - def __repr__(self): - return repr(self.data) + "\n" + repr(self.info) - - def copy(self, copy_data=True): - """Return a copy of the data. - """ - if copy_data: - data = self.data.copy() - else: - data = self.data - return Dataset(data, **self.info) - - -class ProductList(object): - """Class for reading and storing product lists. - """ - def __init__(self, fname): - self.fname = fname - - tree = etree.parse(fname) - self._xml = tree.getroot() - self.prodlist = None - self.attrib = {} - self.vars = {} - self.aliases = {} - self.groups = [] - self.parse() - - def insert_vars(self): - """Variable replacement - """ - for item in self.prodlist.getiterator(): - for key in self.vars: - if key in item.attrib and item.attrib[key] in self.vars[key]: - item.set(key, self.vars[key][item.attrib[key]]) - - def check_groups(self): - """Check area groups. - """ - # create the "rest" group - last_group = [] - for area in self.prodlist: - assigned = False - for group in self.groups: - if len(area.attrib.keys()) == 0: - assigned = True - continue - if area.attrib["id"] in group.data: - assigned = True - break - if not assigned: - last_group.append(area.attrib["id"]) - if last_group: - self.groups.append(Dataset(last_group, id="_rest")) - - # replace area ids with actual xml area items - groups = [] - for group in self.groups: - new_group = Dataset([], **group.info) - for area_id in group.data: - assigned = False - for area in self.prodlist: - if area.attrib["id"] == area_id: - new_group.data.append(area) - assigned = True - break - if not assigned: - LOGGER.warning("Couldn't find area %s in product list", - area_id) - groups.append(new_group) - self.groups = groups - - def parse(self): - """Parse product list XML file. - """ - for item in self._xml: - if item.tag == "product_list": - self.prodlist = item - elif item.tag == "common": - for citem in item: - self.attrib[citem.tag] = citem.text - elif item.tag == "groups": - for group in item: - self.groups.append(Dataset(group.text.split(","), - **group.attrib)) - elif item.tag == "variables": - if item.attrib: - res = [os.environ[env] != val - for env, val in item.attrib.items()] - if any(res): - continue - for var in item: - self.vars.setdefault( - var.tag, {})[var.attrib["id"]] = var.text - elif item.tag == "aliases": - for alias in item: - self.aliases.setdefault( - alias.tag, - {})[alias.attrib["src"]] = alias.attrib['dst'] - self.insert_vars() - self.check_groups() - - -def get_root(fname): - '''Read XML file and return the root tree. - ''' - - tree = etree.parse(fname) - root = tree.getroot() - - return root - - -def parse_xml(tree, also_empty=False): - '''Parse the given XML file to dictionary. - ''' - xml_dict = {} - - # these tags will always be lists, if they are present and non-empty - listify = ['area', 'product', 'valid_satellite', 'invalid_satellite', - 'pattern', 'file_tag', 'directory'] - children = list(tree) - - if len(children) == 0: - try: - xml_dict = tree.text.strip() - except AttributeError: - pass - - for child in children: - new_val = parse_xml(child, also_empty=also_empty) - if len(new_val) == 0: - if also_empty: - xml_dict[child.tag] = '' - continue - else: - continue - - print "child.tag: ", child.tag - - if child.tag in xml_dict: - if not isinstance(xml_dict[child.tag], list): - xml_dict[child.tag] = [xml_dict[child.tag]] - xml_dict[child.tag].append(new_val) - else: - if len(new_val) > 0: - if child.tag in listify: - xml_dict[child.tag] = [new_val] - else: - xml_dict[child.tag] = new_val - - print xml_dict - return xml_dict - - -def get_filepattern_config(fname=None): - '''Retrieves the filepattern configuration file for trollstalker, - and returns the parsed XML as a dictionary. Optional argument - *fname* can be used to specify the file. If *fname* is None, the - systemwide file is read. - ''' - - if fname is None: - fname = os.path.realpath(__file__).split(os.path.sep)[:-2] - fname.append('etc') - fname.append('filepattern_config.xml') - fname = os.path.sep.join(fname) - - return parse_xml(get_root(fname), also_empty=True) From fc723a357b4ba07355d911eaa36d74cd7c16d7de Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Wed, 1 Jun 2016 10:19:30 +0200 Subject: [PATCH 20/27] Bumped version Added IOError on not found xml file --- mpop/satout/mitiff.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 6312fbea..2bbb3a54 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -29,7 +29,7 @@ """mpop mitiff writer interface. """ -__revision__ = 0.1 +__revision__ = 0.2 KELVIN_TO_CELSIUS = -273.15 @@ -42,17 +42,17 @@ logger = logging.getLogger(__name__) -def save(scene, filename, product_name=None): +def save(scene, filename): """Saves the scene as a MITIFF file. This format can be read by DIANA """ - print "inside save ... " - return mitiff_writer(filename, scene, mitiff_product_name=product_name) + logger.debug("inside save ... ") + return mitiff_writer(filename, scene) -def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=None): +def mitiff_writer(filename, root_object, compression=True): """ Write data to MITIFF file. """ - print "Inside mitiff_writer ... " + logger.debug("Inside mitiff_writer ... ") import xml.etree.ElementTree as etree @@ -62,10 +62,9 @@ def mitiff_writer(filename, root_object, compression=True, mitiff_product_name=N if os.path.isfile(fname_): tree = etree.parse(fname_) root = tree.getroot() - - #config = xml_read.parse_xml(xml_read.get_root(fname_)) - #raise ValueError("Could not find a MiTIFF config file ", fname_) - + else: + raise IOError("Could not find %s. Check if this file is available at this location." % (fname_)) + #print "instrument_name: ", root_object.instrument_name product_list = root.findall("product[@instrument_name='%s']" % (root_object.instrument_name)) @@ -272,6 +271,9 @@ def _make_image_description(ro, pc): else: logger.warning("Unknown calib type. Must be Radiance, Reflectance or BT.") + #How to format string by passing the format + #http://stackoverflow.com/questions/1598579/rounding-decimals-with-new-python-format-function + _image_description += ', 8, [ ' for val in range(0,256): #Comma separated list of values @@ -306,7 +308,3 @@ def _make_image_description(ro, pc): save(local_data,"test.mitiff",product_name="AVHRR") print "Complete." - - - - \ No newline at end of file From 970430f64d7986791bbbaf175e5f8de88eff311d Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:40:53 +0000 Subject: [PATCH 21/27] mpop/satout/mitiff.py some more debug info --- mpop/satout/mitiff.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 2bbb3a54..197f4bf0 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -42,7 +42,7 @@ logger = logging.getLogger(__name__) -def save(scene, filename): +def save(scene, filename, **options): """Saves the scene as a MITIFF file. This format can be read by DIANA """ @@ -235,13 +235,13 @@ def _make_image_description(ro, pc): _image_description += ' By: %.6f' % (ro.area.area_extent[3]/1000.) #UR_y _image_description += '\n' - logger.debug("Area extent: ",ro.area.area_extent) + logger.debug("Area extent: {}".format(ro.area.area_extent)) for ch in pc: palette=False #Make calibration. if palette: - logging.debug("Do palette calibration") + logging.debug("Do palette calibration. Not implemented yet.") else: _image_description += 'Table_calibration: ' if 'name-alias' in ch.keys(): From 20ed4717634eee3d12d25269bf747ad02f92f22e Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:41:33 +0000 Subject: [PATCH 22/27] mpop/imageo/geo_image.py more info --- mpop/imageo/geo_image.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/mpop/imageo/geo_image.py b/mpop/imageo/geo_image.py index c1157b03..7b1316da 100644 --- a/mpop/imageo/geo_image.py +++ b/mpop/imageo/geo_image.py @@ -96,6 +96,9 @@ def save(self, filename, compression=6, if fformat.lower() in ('tif', 'tiff'): return self.geotiff_save(filename, compression, tags, gdal_options, blocksize, **kwargs) + elif fformat.lower() in ('mitiff'): + logger.info("Writing to mitiff ....") + try: # Let image.pil_save it ? Image.save(self, filename, compression, fformat=fformat) From 4f6ddd5569afe07cf116898a0afbc71ceec1b7df Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:44:29 +0000 Subject: [PATCH 23/27] mpop/scene.py changed save definition. This could possible break the code. This needs to be confirmed --- mpop/scene.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mpop/scene.py b/mpop/scene.py index 9d993ef1..3fa046fc 100644 --- a/mpop/scene.py +++ b/mpop/scene.py @@ -499,14 +499,16 @@ def load(self, channels=None, load_again=False, area_extent=None, **kwargs): self.channels_to_load = set() - def save(self, filename, to_format="netcdf4", **options): + def save(self, filename, fformat="netcdf4", compression=None, **options): """Saves the current scene into a file of format *to_format*. Supported formats are: - *netcdf4*: NetCDF4 with CF conventions. """ - + to_format = fformat writer = "satout." + to_format + LOG.debug("writer module: %s" %(writer)) + try: writer_module = __import__(writer, globals(), locals(), ["save"]) except ImportError, err: From 5b42afd625c94e5d9049a82cf2465c45ca7ced40 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Fri, 3 Jun 2016 07:46:47 +0000 Subject: [PATCH 24/27] mpop/scene.py added info in mitiff save format --- mpop/scene.py | 1 + 1 file changed, 1 insertion(+) diff --git a/mpop/scene.py b/mpop/scene.py index 3fa046fc..c14127f8 100644 --- a/mpop/scene.py +++ b/mpop/scene.py @@ -504,6 +504,7 @@ def save(self, filename, fformat="netcdf4", compression=None, **options): formats are: - *netcdf4*: NetCDF4 with CF conventions. + - *mitiff*: MITIFF format, TIFF file with image description tag holding the geolocation info. Used by DIANA at Norwegian MET """ to_format = fformat writer = "satout." + to_format From f7d641c089f55e635ef07fb7766d0f9ba0a1fe02 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Thu, 11 Aug 2016 10:04:33 +0000 Subject: [PATCH 25/27] Small fixes and more debug --- mpop/satellites/__init__.py | 5 +++++ mpop/satout/mitiff.py | 19 ++++++++++++++----- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/mpop/satellites/__init__.py b/mpop/satellites/__init__.py index d0d1132b..5a2c2bc2 100644 --- a/mpop/satellites/__init__.py +++ b/mpop/satellites/__init__.py @@ -166,6 +166,7 @@ def create_scene(satname, satnumber, instrument, time_slot, area=None, variant=''): """Create a compound satellite scene. """ + LOG.debug("create_scene") return GenericFactory.create_scene(satname, satnumber, instrument, time_slot, None, area, variant) @@ -182,6 +183,8 @@ def create_scene(satname, satnumber, instrument, time_slot, orbit=None, """Create a compound satellite scene. """ + LOG.debug("create_scene") + return GenericFactory.create_scene(satname, satnumber, instrument, time_slot, orbit, area, variant) @@ -190,6 +193,7 @@ class GenericFactory(object): """Factory for generic satellite scenes. """ + LOG.critical("GenericFactory") @staticmethod def create_scene(satname, satnumber, instrument, time_slot, orbit, @@ -208,6 +212,7 @@ def create_scene(satname, satnumber, instrument, time_slot, orbit, compositer = get_sat_instr_compositer(satellite, instrument) instrument_scene._CompositerClass = compositer + LOG.critical("instrument_scene: {}".format(instrument_scene)) if compositer is not None: # Pass weak ref to compositor to allow garbage collection instrument_scene.image = compositer( diff --git a/mpop/satout/mitiff.py b/mpop/satout/mitiff.py index 197f4bf0..0d4958fb 100644 --- a/mpop/satout/mitiff.py +++ b/mpop/satout/mitiff.py @@ -106,7 +106,9 @@ def mitiff_writer(filename, root_object, compression=True): for ch in product_channels: found_channel = False print ch.attrib['name'] + print root_object.channels for channels in root_object.channels: + print channels if ch.attrib['name'] == channels.name: #Need to scale the data set to mitiff 0-255. 0 is no/missing/bad data. logger.debug("min %f max %f value" % (float(ch.attrib['min-val']),float(ch.attrib['max-val']))) @@ -133,7 +135,12 @@ def mitiff_writer(filename, root_object, compression=True): if not found_channel: logger.debug("Could not find configured channel in read data set. Fill with empty.") - fill_channel = np.zeros(root_object.channels[0].data.shape,dtype=np.uint8) + try: + fill_channel = np.zeros(root_object.channels[0].data.shape,dtype=np.uint8) + except IndexError as ie: + logger.error("Index out out bounds for channels[0].data.shape. Try area instead ... {}".format(ie)) + fill_channel = np.zeros(root_object.area.shape,dtype=np.uint8) + tif.write_image(fill_channel) @@ -216,8 +223,9 @@ def _make_image_description(ro, pc): _image_description += ' +units=km' - _image_description += ' +x_0=%.6f' % (-ro.area.area_extent[0]) - _image_description += ' +y_0=%.6f' % (-ro.area.area_extent[1]) + #Need to use center of lower left pixel. Subtract half a pixel size + _image_description += ' +x_0=%.6f' % (-ro.area.area_extent[0]-ro.area.pixel_size_x/2.) + _image_description += ' +y_0=%.6f' % (-ro.area.area_extent[1]-ro.area.pixel_size_y/2.) _image_description += '\n' _image_description += ' TrueLat: 60N\n' @@ -231,8 +239,9 @@ def _make_image_description(ro, pc): _image_description += ' Ax: %.6f' % (ro.area.pixel_size_x/1000.) _image_description += ' Ay: %.6f' % (ro.area.pixel_size_y/1000.) #But this ads up to upper left corner of upper left pixel. - _image_description += ' Bx: %.6f' % (ro.area.area_extent[0]/1000.) #LL_x - _image_description += ' By: %.6f' % (ro.area.area_extent[3]/1000.) #UR_y + #But need to use the center of the pixel. Therefor use the center of the upper left pixel. + _image_description += ' Bx: %.6f' % (ro.area.area_extent[0]/1000. + ro.area.pixel_size_x/1000./2.) #LL_x + _image_description += ' By: %.6f' % (ro.area.area_extent[3]/1000. - ro.area.pixel_size_y/1000./2.) #UR_y _image_description += '\n' logger.debug("Area extent: {}".format(ro.area.area_extent)) From a61ceae92b5c1c0d68ae857589490f21e4c45619 Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Mon, 14 Nov 2016 12:39:44 +0000 Subject: [PATCH 26/27] Removed debugging output --- mpop/satellites/__init__.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/mpop/satellites/__init__.py b/mpop/satellites/__init__.py index 5a2c2bc2..db7d7690 100644 --- a/mpop/satellites/__init__.py +++ b/mpop/satellites/__init__.py @@ -166,7 +166,6 @@ def create_scene(satname, satnumber, instrument, time_slot, area=None, variant=''): """Create a compound satellite scene. """ - LOG.debug("create_scene") return GenericFactory.create_scene(satname, satnumber, instrument, time_slot, None, area, variant) @@ -183,8 +182,6 @@ def create_scene(satname, satnumber, instrument, time_slot, orbit=None, """Create a compound satellite scene. """ - LOG.debug("create_scene") - return GenericFactory.create_scene(satname, satnumber, instrument, time_slot, orbit, area, variant) @@ -193,8 +190,6 @@ class GenericFactory(object): """Factory for generic satellite scenes. """ - LOG.critical("GenericFactory") - @staticmethod def create_scene(satname, satnumber, instrument, time_slot, orbit, area=None, variant=''): @@ -212,7 +207,6 @@ def create_scene(satname, satnumber, instrument, time_slot, orbit, compositer = get_sat_instr_compositer(satellite, instrument) instrument_scene._CompositerClass = compositer - LOG.critical("instrument_scene: {}".format(instrument_scene)) if compositer is not None: # Pass weak ref to compositor to allow garbage collection instrument_scene.image = compositer( From b3df972f3f4a976cc89f34a4eca423ac2a5fa86a Mon Sep 17 00:00:00 2001 From: Trygve Aspenes Date: Tue, 15 Nov 2016 12:23:17 +0100 Subject: [PATCH 27/27] Try to fix conflict --- mpop/imageo/geo_image.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mpop/imageo/geo_image.py b/mpop/imageo/geo_image.py index 7b1316da..df6cfa09 100644 --- a/mpop/imageo/geo_image.py +++ b/mpop/imageo/geo_image.py @@ -95,7 +95,8 @@ def save(self, filename, compression=6, if fformat.lower() in ('tif', 'tiff'): return self.geotiff_save(filename, compression, tags, - gdal_options, blocksize, **kwargs) + gdal_options, blocksize, + **kwargs) elif fformat.lower() in ('mitiff'): logger.info("Writing to mitiff ....")