# -*- coding: utf-8 -*-
"""
ninjotiff.py
Created on Mon Apr 15 13:41:55 2013
A big amount of the tiff writer are (PFE) from
https://github.com/davidh-ssec/polar2grid by David Hoese
License:
Copyright (C) 2013 Space Science and Engineering Center (SSEC),
University of Wisconsin-Madison.
Lars Ørum Rasmussen, DMI.
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 .
Original scripts and automation included as part of this package are
distributed under the GNU GENERAL PUBLIC LICENSE agreement version 3.
Binary executable files included as part of this software package are
copyrighted and licensed by their respective organizations, and
distributed consistent with their licensing terms.
"""
import os
import copy
import logging
import calendar
from datetime import datetime
import numpy as np
import mpop.imageo.formats.libtiff as libtiff
from mpop.imageo.formats.libtiff import TIFF, TIFFFieldInfo, TIFFDataType, FIELD_CUSTOM
log = logging.getLogger(__name__)
#-------------------------------------------------------------------------------
#
# Ninjo tiff tags from DWD
#
#-------------------------------------------------------------------------------
# Geotiff tags
GTF_ModelPixelScale = 33550
GTF_ModelTiepoint = 33922
NTD_Magic = 40000
NTD_SatelliteNameID = 40001
NTD_DateID = 40002
NTD_CreationDateID = 40003
NTD_ChannelID = 40004
NTD_HeaderVersion = 40005
NTD_FileName = 40006
NTD_DataType = 40007
NTD_SatelliteNumber = 40008
NTD_ColorDepth = 40009
NTD_DataSource = 40010
NTD_XMinimum = 40011
NTD_XMaximum = 40012
NTD_YMinimum = 40013
NTD_YMaximum = 40014
NTD_Projection = 40015
NTD_MeridianWest = 40016
NTD_MeridianEast = 40017
NTD_EarthRadiusLarge = 40018
NTD_EarthRadiusSmall = 40019
NTD_GeodeticDate = 40020
NTD_ReferenceLatitude1 = 40021
NTD_ReferenceLatitude2 = 40022
NTD_CentralMeridian = 40023
NTD_PhysicValue = 40024
NTD_PhysicUnit = 40025
NTD_MinGrayValue = 40026
NTD_MaxGrayValue = 40027
NTD_Gradient = 40028
NTD_AxisIntercept = 40029
NTD_ColorTable = 40030
NTD_Description = 40031
NTD_OverflightDirection = 40032
NTD_GeoLatitude = 40033
NTD_GeoLongitude = 40034
NTD_Altitude = 40035
NTD_AOSAsimuth = 40036
NTD_LOSAsimuth = 40037
NTD_MaxElevation = 40038
NTD_OverflightTime = 40039
NTD_IsBlackLineCorrection = 40040
NTD_IsAtmosphereCorrected = 40041
NTD_IsCalibrated = 40042
NTD_IsNormalized = 40043
NTD_OriginalHeader = 40044
NTD_IsValueTableAvailable = 40045
NTD_ValueTableStringField = 40046
NTD_ValueTableFloatField = 40047
NTD_TransparentPixel = 50000
#
# model_pixel_scale_tag_count ? ...
# Sometimes DWD product defines an array of length 2 (instead of 3 (as in geotiff)).
#
MODEL_PIXEL_SCALE_COUNT = int(os.environ.get("GEOTIFF_MODEL_PIXEL_SCALE_COUNT", 3))
ninjo_tags_dict = {
# Geotiff tags
GTF_ModelPixelScale:
TIFFFieldInfo(GTF_ModelPixelScale, MODEL_PIXEL_SCALE_COUNT,
MODEL_PIXEL_SCALE_COUNT, TIFFDataType.TIFF_DOUBLE,
FIELD_CUSTOM, True, False, "ModelPixelScale" ),
GTF_ModelTiepoint:
TIFFFieldInfo(GTF_ModelTiepoint, 6, 6, TIFFDataType.TIFF_DOUBLE,
FIELD_CUSTOM, True, False, "ModelTiePoint" ),
# DWD tags
NTD_Magic:
TIFFFieldInfo(NTD_Magic, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "Magic" ),
NTD_SatelliteNameID:
TIFFFieldInfo(NTD_SatelliteNameID, 1, 1, TIFFDataType.TIFF_LONG,
FIELD_CUSTOM, True, False, "SatelliteNameID" ),
NTD_DateID:
TIFFFieldInfo(NTD_DateID, 1, 1, TIFFDataType.TIFF_LONG,
FIELD_CUSTOM, True, False, "DateID" ),
NTD_CreationDateID:
TIFFFieldInfo(NTD_CreationDateID, 1, 1, TIFFDataType.TIFF_LONG,
FIELD_CUSTOM, True, False, "CreationDateID" ),
NTD_ChannelID:
TIFFFieldInfo(NTD_ChannelID, 1, 1, TIFFDataType.TIFF_LONG,
FIELD_CUSTOM, True, False, "ChannelID" ),
NTD_HeaderVersion:
TIFFFieldInfo(NTD_HeaderVersion, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "HeaderVersion" ),
NTD_FileName:
TIFFFieldInfo(NTD_FileName, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "FileName" ),
NTD_DataType:
TIFFFieldInfo(NTD_DataType, 5, 5, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "DataType" ), # 4 chars + NUL character
NTD_SatelliteNumber:
TIFFFieldInfo(NTD_SatelliteNumber, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "SatelliteNumber" ),
NTD_ColorDepth:
TIFFFieldInfo(NTD_ColorDepth, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "ColorDepth" ),
NTD_DataSource:
TIFFFieldInfo(NTD_DataSource, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "DataSource" ),
NTD_XMinimum:
TIFFFieldInfo(NTD_XMinimum, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "XMinimum" ),
NTD_XMaximum:
TIFFFieldInfo(NTD_XMaximum, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "XMaximum" ),
NTD_YMinimum:
TIFFFieldInfo(NTD_YMinimum, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "YMinimum" ),
NTD_YMaximum:
TIFFFieldInfo(NTD_YMaximum, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "YMaximum" ),
NTD_Projection:
TIFFFieldInfo(NTD_Projection, 5, 5, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "Projection" ), # 4 chars + NUL character
NTD_MeridianWest:
TIFFFieldInfo(NTD_MeridianWest, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "MeridianWest" ),
NTD_MeridianEast:
TIFFFieldInfo(NTD_MeridianEast, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "MeridianEast" ),
NTD_EarthRadiusLarge:
TIFFFieldInfo(NTD_EarthRadiusLarge, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "EarthRadiusLarge" ),
NTD_EarthRadiusSmall:
TIFFFieldInfo(NTD_EarthRadiusSmall, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "EarthRadiusSmall" ),
NTD_GeodeticDate:
TIFFFieldInfo(NTD_GeodeticDate, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "GeodeticDate" ), # Max 20
NTD_ReferenceLatitude1:
TIFFFieldInfo(NTD_ReferenceLatitude1, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "ReferenceLatitude1" ),
NTD_ReferenceLatitude2:
TIFFFieldInfo(NTD_ReferenceLatitude2, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "ReferenceLatitude2" ),
NTD_CentralMeridian:
TIFFFieldInfo(NTD_CentralMeridian, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "CentralMeridian" ),
NTD_PhysicValue:
TIFFFieldInfo(NTD_PhysicValue, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "PhysicValue" ), # Max 10
NTD_PhysicUnit:
TIFFFieldInfo(NTD_PhysicUnit, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "PhysicUnit" ), # Max 10
NTD_MinGrayValue:
TIFFFieldInfo(NTD_MinGrayValue, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "MinGrayValue" ),
NTD_MaxGrayValue:
TIFFFieldInfo(NTD_MaxGrayValue, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "MaxGrayValue" ),
NTD_Gradient:
TIFFFieldInfo(NTD_Gradient, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "Gradient" ),
NTD_AxisIntercept:
TIFFFieldInfo(NTD_AxisIntercept, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "AxisIntercept" ),
NTD_ColorTable:
TIFFFieldInfo(NTD_ColorTable, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "ColorTable" ),
NTD_Description:
TIFFFieldInfo(NTD_Description, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "Description" ),
NTD_OverflightDirection:
TIFFFieldInfo(NTD_OverflightDirection, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "OverflightDirection" ),
NTD_GeoLatitude:
TIFFFieldInfo(NTD_GeoLatitude, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "GeoLatitude" ),
NTD_GeoLongitude:
TIFFFieldInfo(NTD_GeoLongitude, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "GeoLongitude" ),
NTD_Altitude:
TIFFFieldInfo(NTD_Altitude, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "Altitude" ),
NTD_AOSAsimuth:
TIFFFieldInfo(NTD_AOSAsimuth, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "AOSAsimuth" ),
NTD_LOSAsimuth:
TIFFFieldInfo(NTD_LOSAsimuth, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "LOSAsimuth" ),
NTD_MaxElevation:
TIFFFieldInfo(NTD_MaxElevation, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "MaxElevation" ),
NTD_OverflightTime:
TIFFFieldInfo(NTD_OverflightTime, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "OverflightTime" ),
NTD_IsBlackLineCorrection:
TIFFFieldInfo(NTD_IsBlackLineCorrection, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "IsBlackLineCorrection" ),
NTD_IsAtmosphereCorrected:
TIFFFieldInfo(NTD_IsAtmosphereCorrected, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "IsAtmosphereCorrected" ),
NTD_IsCalibrated:
TIFFFieldInfo(NTD_IsCalibrated, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "IsCalibrated" ),
NTD_IsNormalized:
TIFFFieldInfo(NTD_IsNormalized, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "IsNormalized" ),
NTD_OriginalHeader:
TIFFFieldInfo(NTD_OriginalHeader, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "OriginalHeader" ),
NTD_IsValueTableAvailable:
TIFFFieldInfo(NTD_IsValueTableAvailable, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "IsValueTableAvailable" ),
NTD_ValueTableStringField:
TIFFFieldInfo(NTD_ValueTableStringField, -1, -1, TIFFDataType.TIFF_ASCII,
FIELD_CUSTOM, True, False, "ValueTableStringField" ),
NTD_ValueTableFloatField:
TIFFFieldInfo(NTD_ValueTableFloatField, 1, 1, TIFFDataType.TIFF_FLOAT,
FIELD_CUSTOM, True, False, "ValueTableFloatField" ),
NTD_TransparentPixel:
TIFFFieldInfo(NTD_TransparentPixel, 1, 1, TIFFDataType.TIFF_SLONG,
FIELD_CUSTOM, True, False, "TransparentPixel" ),
}
# Add Ninjo tags to the libtiff library
_ninjo_tags_extender = libtiff.add_tags(ninjo_tags_dict.values())
ninjo_tags = sorted(ninjo_tags_dict.keys())
#-------------------------------------------------------------------------------
#
# Read Ninjo products config file.
#
#-------------------------------------------------------------------------------
def get_product_config(product_name, force_read=False):
"""Read Ninjo configuration entry for a given product name.
:Parameters:
product_name : str
Name of Ninjo product.
:Arguments:
force_read : Boolean
Force re-reading config file.
**Notes**:
* It will look for a *ninjotiff_products.cfg* in MPOP's
configuration directory defined by *PPP_CONFIG_DIR*.
* As an example, see *ninjotiff_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()
#print "Reading Ninjo config file: '%s'" % filename
log.info("Reading Ninjo 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)
products[sec] = prd
self._products = products
@staticmethod
def _find_a_config_file():
name_ = 'ninjotiff_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 Ninjo tiff config file")
#-------------------------------------------------------------------------------
#
# Read tiff file.
#
#-------------------------------------------------------------------------------
class _TIFF(object):
""" Just an context wrapper around an libtiff.TIFF instance.
"""
def __init__(self, filename, mode='r'):
"""Open a tiff file.
see: libtiff.TIFF.open()
"""
self.tiff = TIFF.open(filename, mode)
self.tiff.ninjo_tags_dict = ninjo_tags_dict
self.tiff.ninjo_tags = ninjo_tags
def __enter__(self):
return self.tiff
def __exit__(self, type_, value, traceback):
self.tiff.close()
def _read_directories(self):
"""Iterate over directories in a tiff file.
:Parameters:
self : libtiff.TIFF
A TIFF instance.
:Returns:
tiff_directory : Tiff object
A Tiff directory instance.
"""
yield self
while not self.LastDirectory():
self.ReadDirectory()
yield self
self.SetDirectory(0)
def info(filename):
"""Read metadata from Tiff file.
:Parameters:
filename : str
Name of Tiff file.
:Returns:
iterator : a Python generator iterator
A "list" of tiff metadata.
**Usage**::
for inf in info(filename):
print inf, '\n'
"""
with _TIFF(filename) as self:
for d in _read_directories(self):
l = []
for item in d.info().split('\n'):
k, v = item.split(':', 1)
if (k.endswith('OffSets') or
k.endswith('ByteCounts') or
k == 'FileName' or
k == 'DataType'):
continue
l.append(item)
for tag in d.ninjo_tags:
value = d.GetField(tag)
name = d.ninjo_tags_dict[tag].field_name
if value is None:
continue
l.append('%s: %s' % (name, str(value)))
yield '\n'.join(l)
def image_data(filename):
"""Read image data from Tiff file.
:Parameters:
filename : str
Name of Tiff file.
**Usage**::
for img in image_data(filename):
print img
"""
with _TIFF(filename) as self:
for d in _read_directories(self):
yield d.read_tiles()
def colortable(filename):
"""Read colortables from Tiff file.
:Parameters:
filename : str
Name of Tiff file.
**Usage**::
for clt in colortable(filename):
print clt
"""
with _TIFF(filename) as self:
return self.GetField('ColorMap')
#-------------------------------------------------------------------------------
#
# Write Ninjo Products
#
#-------------------------------------------------------------------------------
def _get_physic_value(physic_unit):
# return Ninjo's physics unit and value.
if physic_unit.upper() in ('K', 'KELVIN'):
return 'Kelvin', 'T'
elif physic_unit.upper() in ('C', 'CELSIUS'):
return 'Celsius', 'T'
elif physic_unit == '%':
return physic_unit, 'Reflectance'
elif physic_unit.upper() in ('MW M-2 SR-1 (CM-1)-1',):
return physic_unit, 'Radiance'
else:
return physic_unit, 'Unknown'
def _get_projection_name(area_def):
# return Ninjo's projection name.
proj_name = area_def.proj_dict['proj']
if proj_name in ('eqc',):
return 'PLAT'
elif proj_name in ('stere',):
lat_0 = area_def.proj_dict['lat_0']
if lat_0 < 0:
return 'SPOL'
else:
return 'NPOL'
return None
def _finalize(geo_image):
"""Finalize a mpop GeoImage for Ninjo. Specialy take care of phycical scale
and offset.
:Parameters:
geo_image : mpop.imageo.geo_image.GeoImage
See MPOP's documentation.
:Returns:
image : numpy.array
Final image.
scale : float
Scale for transform pixel value to physical value.
offset : float
Offset for transform pixel value to physical value.
fill_value : int
Value for used masked out pixels.
**Notes**:
physic_val = image*scale + offset
"""
if geo_image.mode == 'L':
# PFE: mpop.satout.cfscene
dtype = np.uint8
data = geo_image.channels[0]
fill_value = geo_image.fill_value or 0
if np.ma.count_masked(data) == data.size:
# All data is masked
data = np.ones(data.shape, dtype=dtype) * fill_value
scale = 1
offset = 0
else:
chn_max = data.max()
chn_min = data.min()
scale = ((chn_max - chn_min) /
(2**np.iinfo(dtype).bits - 1.0))
# Handle the case where all data has the same value.
scale = scale or 1
offset = chn_min
mask = data.mask
data = ((data.data - offset) / scale).astype(dtype)
data[mask] = fill_value
return data, scale, offset, fill_value
elif geo_image.mode == 'RGB':
channels, fill_value = geo_image._finalize()
fill_value = fill_value or (0, 0, 0)
data = np.dstack((channels[0].filled(fill_value[0]),
channels[1].filled(fill_value[1]),
channels[2].filled(fill_value[2])))
return data, 1.0, 0.0, fill_value[0]
else:
raise ValueError("Don't known how til handle image mode '%s'" %
str(geo_image.mode))
def save(geo_image, filename, ninjo_product_name=None, **kwargs):
"""MPOP's interface to Ninjo TIFF writer.
:Parameters:
geo_image : mpop.imageo.geo_image.GeoImage
See MPOP's documentation.
filename : str
The name of the TIFF file to be created
:Keywords:
ninjo_product_name : str
Optional index to Ninjo configuration file.
kwargs : dict
See _write
"""
data, scale, offset, fill_value = _finalize(geo_image)
area_def = geo_image.area
time_slot = geo_image.time_slot
# Some Ninjo tiff names
kwargs['image_dt'] = time_slot
kwargs['transparent_pix'] = fill_value
kwargs['gradient'] = scale
kwargs['axis_intercept'] = offset
kwargs['is_calibrated'] = True
write(data, filename, area_def, ninjo_product_name, **kwargs)
def write(image_data, output_fn, area_def, product_name=None, **kwargs):
"""Generic Ninjo TIFF writer.
If 'prodcut_name' is given, it will load corresponding Ninjo tiff metadata
from '${PPP_CONFIG_DIR}/ninjotiff.cfg'. Else, all Ninjo tiff metadata should
be passed by '**kwargs'. A mixture is allowed, where passed arguments
overwrite config file.
:Parameters:
image_data : 2D numpy array
Satellite image data to be put into the NinJo compatible tiff
output_fn : str
The name of the TIFF file to be created
area_def: pyresample.geometry.AreaDefinition
Defintion of area
product_name : str
Optional index to Ninjo configuration file.
:Keywords:
kwargs : dict
See _write
"""
upper_left = area_def.get_lonlat(0, 0)
lower_right = area_def.get_lonlat(area_def.shape[0], area_def.shape[1])
scale = abs(lower_right[0] - upper_left[0])/area_def.shape[1],\
abs(upper_left[1] - lower_right[1])/area_def.shape[0]
if len(image_data.shape) == 3:
shape = (area_def.y_size, area_def.x_size, 3)
write_rgb = True
log.info("Will generate RGB product '%s'" % product_name)
else:
shape = (area_def.y_size, area_def.x_size)
write_rgb = False
log.info("Will generate product '%s'" % product_name)
if image_data.shape != shape:
raise ValueError, "Raster shape %s does not correspond to expected shape %s" % (
str(image_data.shape), str(shape))
# Ninjo's physical units and value.
# If just a physical unit (e.g. 'C') is passed, it will then be
# translated into Ninjo's unit and value (e.q 'CELCIUS' and 'T').
physic_unit = kwargs.get('physic_unit', None)
if physic_unit and not kwargs.get('physic_value', None):
kwargs['physic_unit'], kwargs['physic_value'] = \
_get_physic_value(physic_unit)
# Ninjo's projection name.
kwargs['projection'] = kwargs.pop('projection', None) or \
_get_projection_name(area_def) or \
area_def.proj_id.split('_')[-1]
if product_name:
options = get_product_config(product_name)
else:
options = {}
options['meridian_west'] = upper_left[0]
options['meridian_east'] = lower_right[0]
options['pixel_xres'] = scale[0]
options['pixel_yres'] = scale[1]
options['origin_lon'] = upper_left[0]
options['origin_lat'] = upper_left[1]
options['min_gray_val'] = image_data.min()
options['max_gray_val'] = image_data.max()
options.update(kwargs) # Update/overwrite with passed arguments
_write(image_data, output_fn, write_rgb=write_rgb, **options)
#-------------------------------------------------------------------------------
#
# Write tiff file.
#
#-------------------------------------------------------------------------------
def _write(image_data, output_fn, write_rgb=False, **kwargs):
"""Proudly Found Elsewhere (PFE) https://github.com/davidh-ssec/polar2grid
by David Hoese.
Create a NinJo compatible TIFF file with the tags used
by the DWD's version of NinJo. Also stores the image as tiles on disk
and creates a multi-resolution/pyramid/overview set of images
(deresolution: 2,4,8,16).
:Parameters:
image_data : 2D or 3D numpy array
Satellite image data to be put into the NinJo compatible tiff
An 3D array (HxWx3) is expected for a RGB image.
filename : str
The name of the TIFF file to be created
:Keywords:
cmap : tuple/list of 3 lists of uint16's
Individual RGB arrays describing the color value for the
corresponding data value. For example, image data with a data
type of unsigned 8-bit integers have 256 possible values (0-255).
So each list in cmap will have 256 values ranging from 0 to
65535 (2**16 - 1). (default linear B&W colormap)
sat_id : int
DWD NinJo Satellite ID number
chan_id : int
DWD NinJo Satellite Channel ID number
data_source : str
String describing where the data came from (SSEC, EUMCAST)
tile_width : int
Width of tiles on disk (default 512)
tile_length : int
Length of tiles on disk (default 512)
data_cat : str
NinJo specific data category
- data_cat[0] = P (polar) or G (geostat)
- data_cat[1] = O (original) or P (product)
- data_cat[2:4] = RN or RB or RA or RN or AN (Raster, Bufr, ASCII, NIL)
Example: 'PORN' or 'GORN' or 'GPRN' or 'PPRN'
pixel_xres : float
Nadir view pixel resolution in degrees longitude
pixel_yres : float
Nadir view pixel resolution in degrees latitude
origin_lat : float
Top left corner latitude
origin_lon : float
Top left corner longitude
image_dt : datetime object
Python datetime object describing the date and time of the image
data provided in UTC
projection : str
NinJo compatible projection name (NPOL,PLAT,etc.)
meridian_west : float
Western image border (default 0.0)
meridian_east : float
Eastern image border (default 0.0)
radius_a : float
Large/equatorial radius of the earth (default )
radius_b : float
Small/polar radius of the earth (default )
ref_lat1 : float
Reference latitude 1 (default )
ref_lat2 : float
Reference latitude 2 (default )
central_meridian : float
Central Meridian (default )
physic_value : str
Physical value type. Examples:
- Temperature = 'T'
- Albedo = 'ALBEDO'
physic_unit : str
Physical value units. Examples:
- 'CELSIUS'
- '%'
min_gray_val : int
Minimum gray value (default 0)
max_gray_val : int
Maximum gray value (default 255)
gradient : float
Gradient/Slope
axis_intercept : float
Axis Intercept
altitude : float
Altitude of the data provided (default 0.0)
is_atmo_corrected : bool
Is the data atmosphere corrected? (True/1 for yes) (default False/0)
is_calibrated : bool
Is the data calibrated? (True/1 for yes) (default False/0)
is_normalized : bool
Is the data normalized (True/1 for yes) (default False/0)
description : str
Description string to be placed in the output TIFF (optional)
transparent_pix : int
Transparent pixel value (default -1)
:Raises:
KeyError :
if required keyword is not provided
"""
def _raise_value_error(text):
log.error(text)
raise ValueError(text)
def _default_colormap(reverse=False):
# Basic B&W colormap
if reverse:
return [[ x*256 for x in range(255, -1, -1) ]]*3
return [[ x*256 for x in range(256) ]]*3
def _eval_or_none(key, eval_func):
try:
return eval_func(kwargs[key])
except KeyError:
return None
log.info("Creating output file '%s'" % (output_fn,))
tiff = TIFF.open(output_fn, "w")
# Extract keyword arguments
cmap = kwargs.pop("cmap", None)
sat_id = int(kwargs.pop("sat_id"))
chan_id = int(kwargs.pop("chan_id"))
data_source = str(kwargs.pop("data_source"))
tile_width = int(kwargs.pop("tile_width", 512))
tile_length = int(kwargs.pop("tile_length", 512))
data_cat = str(kwargs.pop("data_cat"))
pixel_xres = float(kwargs.pop("pixel_xres"))
pixel_yres = float(kwargs.pop("pixel_yres"))
origin_lat = float(kwargs.pop("origin_lat"))
origin_lon = float(kwargs.pop("origin_lon"))
image_dt = kwargs.pop("image_dt")
projection = str(kwargs.pop("projection"))
meridian_west = float(kwargs.pop("meridian_west", 0.0))
meridian_east = float(kwargs.pop("meridian_east", 0.0))
radius_a = _eval_or_none("radius_a", float)
radius_b = _eval_or_none("radius_b", float)
ref_lat1 = _eval_or_none("ref_lat1", float)
ref_lat2 = _eval_or_none("ref_lat2", float)
central_meridian = _eval_or_none("central_meridian", float)
min_gray_val = int(kwargs.pop("min_gray_val", 0))
max_gray_val = int(kwargs.pop("max_gray_val", 255))
altitude = float(kwargs.pop("altitude", 0.0))
is_blac_corrected = int(bool(kwargs.pop("is_blac_corrected", 0)))
is_atmo_corrected = int(bool(kwargs.pop("is_atmo_corrected", 0)))
is_calibrated = int(bool(kwargs.pop("is_calibrated", 0)))
is_normalized = int(bool(kwargs.pop("is_normalized", 0)))
description = _eval_or_none("description", str)
physic_value = str(kwargs.pop("physic_value", 'None'))
physic_unit = str(kwargs.pop("physic_unit", 'None'))
gradient = float(kwargs.pop("gradient", 1.0))
axis_intercept = float(kwargs.pop("axis_intercept", 0.0))
transparent_pix = int(kwargs.pop("transparent_pix", -1))
# Keyword checks / verification
if not cmap:
if physic_value == 'T':
reverse = True
else:
reverse = False
cmap = _default_colormap(reverse)
if len(cmap) != 3:
_raise_value_error("Colormap (cmap) must be a list of 3 lists (RGB), not %d" %
len(cmap))
if len(data_cat) != 4:
_raise_value_error("NinJo data type must be 4 characters")
if data_cat[0] not in ["P", "G"]:
_raise_value_error("NinJo data type's first character must be 'P' or 'G' not '%s'" %
data_cat[0])
if data_cat[1] not in ["O", "P"]:
_raise_value_error("NinJo data type's second character must be 'O' or 'P' not '%s'" %
data_cat[1])
if data_cat[2:4] not in ["RN","RB","RA","BN","AN"]:
_raise_value_error("NinJo data type's last 2 characters must be one of %s not '%s'" %
("['RN','RB','RA','BN','AN']", data_cat[2:4]))
if description is not None and len(description) >= 1000:
log.error("NinJo description must be less than 1000 characters")
raise ValueError("NinJo description must be less than 1000 characters")
file_dt = datetime.utcnow()
file_epoch = calendar.timegm(file_dt.timetuple())
image_epoch = calendar.timegm(image_dt.timetuple())
def _write_oneres(image_data, pixel_xres, pixel_yres):
log.info("Writing tags and data for a resolution %dx%d" % image_data.shape[:2])
# Write Tag Data
# Built ins
tiff.SetField("ImageWidth", image_data.shape[1])
tiff.SetField("ImageLength", image_data.shape[0])
tiff.SetField("BitsPerSample", 8)
tiff.SetField("Compression", libtiff.COMPRESSION_DEFLATE)
if write_rgb:
tiff.SetField("Photometric", libtiff.PHOTOMETRIC_RGB)
tiff.SetField("SamplesPerPixel", 3)
else:
tiff.SetField("Photometric", libtiff.PHOTOMETRIC_PALETTE)
tiff.SetField("SamplesPerPixel", 1)
tiff.SetField("ColorMap", cmap)
tiff.SetField("Orientation", libtiff.ORIENTATION_TOPLEFT)
tiff.SetField("SMinSampleValue", 0)
tiff.SetField("SMaxsampleValue", 255)
tiff.SetField("PlanarConfig", libtiff.PLANARCONFIG_CONTIG)
tiff.SetField("TileWidth", tile_width)
tiff.SetField("TileLength", tile_length)
tiff.SetField("SampleFormat", libtiff.SAMPLEFORMAT_UINT)
# NinJo specific tags
if description is not None:
tiff.SetField("Description", description)
if MODEL_PIXEL_SCALE_COUNT == 3:
tiff.SetField("ModelPixelScale", [pixel_xres, pixel_yres, 0.0])
else:
tiff.SetField("ModelPixelScale", [pixel_xres, pixel_yres])
tiff.SetField("ModelTiePoint", [0.0, 0.0, 0.0, origin_lon, origin_lat, 0.0])
tiff.SetField("Magic", "NINJO")
tiff.SetField("SatelliteNameID", sat_id)
tiff.SetField("DateID", image_epoch)
tiff.SetField("CreationDateID", file_epoch)
tiff.SetField("ChannelID", chan_id)
tiff.SetField("HeaderVersion", 2)
tiff.SetField("FileName", output_fn)
tiff.SetField("DataType", data_cat)
tiff.SetField("SatelliteNumber", "\x00") # Hardcoded to 0
if write_rgb:
tiff.SetField("ColorDepth", 24)
elif cmap:
tiff.SetField("ColorDepth", 16)
else:
tiff.SetField("ColorDepth", 8)
tiff.SetField("DataSource", data_source)
tiff.SetField("XMinimum", 1)
tiff.SetField("XMaximum", image_data.shape[1])
tiff.SetField("YMinimum", 1)
tiff.SetField("YMaximum", image_data.shape[0])
tiff.SetField("Projection", projection)
tiff.SetField("MeridianWest", meridian_west)
tiff.SetField("MeridianEast", meridian_east)
if radius_a is not None:
tiff.SetField("EarthRadiusLarge", float(radius_a))
if radius_b is not None:
tiff.SetField("EarthRadiusSmall", float(radius_b))
#tiff.SetField("GeodeticDate", "\x00") # ---?
if ref_lat1 is not None:
tiff.SetField("ReferenceLatitude1", ref_lat1)
if ref_lat2 is not None:
tiff.SetField("ReferenceLatitude2", ref_lat2)
if central_meridian is not None:
tiff.SetField("CentralMeridian", central_meridian)
tiff.SetField("PhysicValue", physic_value)
tiff.SetField("PhysicUnit", physic_unit)
tiff.SetField("MinGrayValue", min_gray_val)
tiff.SetField("MaxGrayValue", max_gray_val)
tiff.SetField("Gradient", gradient)
tiff.SetField("AxisIntercept", axis_intercept)
tiff.SetField("Altitude", altitude)
tiff.SetField("IsBlackLineCorrection", is_blac_corrected)
tiff.SetField("IsAtmosphereCorrected", is_atmo_corrected)
tiff.SetField("IsCalibrated", is_calibrated)
tiff.SetField("IsNormalized", is_normalized)
tiff.SetField("TransparentPixel", transparent_pix)
# Write Base Data Image
tiff.write_tiles(image_data)
tiff.WriteDirectory()
# Write multi-resolution overviews (or not)
tiff.SetDirectory(0)
_write_oneres(image_data, pixel_xres, pixel_yres)
for index, scale in enumerate((2, 4, 8, 16)):
shape = (image_data.shape[0]/scale,
image_data.shape[1]/scale)
if shape[0] > tile_width and shape[1] > tile_length:
tiff.SetDirectory(index + 1)
_write_oneres(image_data[::scale,::scale], pixel_xres*scale, pixel_yres*scale)
tiff.close()
log.info("Successfully created a NinJo tiff file: '%s'" % (output_fn,))
if __name__ == '__main__':
import sys
try:
filename = sys.argv[1]
except IndexError:
print >> sys.stderr, "usage: python ninjotiff.py "
sys.exit(2)
for inf in info(filename):
print inf, '\n'