Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature mitiff palette #947

Merged
merged 8 commits into from
Nov 4, 2019
Merged

Feature mitiff palette #947

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
d9d3dad
Add writing of mitiff palette products. Adjustments to projstring han…
Oct 9, 2019
681b665
Dont add aditional information in the writer
Oct 11, 2019
c62fca5
Merge branch 'master' of github.com:pytroll/satpy into feature-mitiff…
Oct 31, 2019
13212e1
add test for palette image
Nov 1, 2019
57210c8
Need to handle when x_0=0 and x_y=0 is set when datasets is a list.
Nov 1, 2019
b0a1e4b
stickler
Nov 1, 2019
d26a1b3
codefactor
Nov 1, 2019
a2375bc
more codefactor
Nov 1, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 87 additions & 0 deletions satpy/tests/writer_tests/test_mitiff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -691,6 +691,93 @@ def test_convert_proj4_string(self):
proj4_string = w._add_proj4_string(ds1, ds1)
self.assertEqual(proj4_string, check['proj4'])

def test_save_dataset_palette(self):
"""Test writer operation as palette."""
import os
import numpy as np
from libtiff import TIFF
from satpy.writers.mitiff import MITIFFWriter

expected = np.full((100, 200), 0)

exp_c = ([0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[2, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

color_map = [[0, 3], [1, 4], [2, 5]]
pal_desc = ['test', 'test2']
unit = "Test"

dataset = self._get_test_one_dataset()
palette = {'palette': True,
'palette_color_map': color_map,
'palette_description': pal_desc,
'palette_unit': unit,
'palette_channel_name': dataset.attrs['name']}

w = MITIFFWriter(base_dir=self.base_dir)
w.save_dataset(dataset, **palette)
filename = "{:s}_{:%Y%m%d_%H%M%S}.mitiff".format(dataset.attrs['name'],
dataset.attrs['start_time'])
tif = TIFF.open(os.path.join(self.base_dir, filename))
# Need to check PHOTOMETRIC is 3, ie palette
self.assertEqual(tif.GetField('PHOTOMETRIC'), 3)
colormap = tif.GetField('COLORMAP')
# Check the colormap of the palette image
self.assertEqual(colormap, exp_c)
IMAGEDESCRIPTION = 270
imgdesc = (tif.GetField(IMAGEDESCRIPTION)).decode('utf-8').split('\n')
found_color_info = False
unit_name_found = False
name_length_found = False
name_length = 0
names = []
unit_name = None
for key in imgdesc:
if name_length_found and name_length > len(names):
names.append(key)
continue
elif unit_name_found:
name_length = int(key)
name_length_found = True
unit_name_found = False
elif found_color_info:
unit_name = key
unit_name_found = True
found_color_info = False
elif 'COLOR INFO:' in key:
found_color_info = True
# Check the name of the palette description
self.assertEqual(name_length, 2)
# Check the name and unit name of the palette
self.assertEqual(unit_name, ' Test')
# Check the palette description of the palette
self.assertEqual(names, [' test', ' test2'])
for image in tif.iter_images():
np.testing.assert_allclose(image, expected, atol=1.e-6, rtol=0)


def suite():
"""The test suite for this writer's tests.
Expand Down
107 changes: 88 additions & 19 deletions satpy/writers/mitiff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2018 Satpy developers
# Copyright (c) 2018, 2019 Satpy developers
#
# This file is part of satpy.
#
Expand All @@ -20,7 +20,6 @@
"""

import logging

import numpy as np

from satpy.writers import ImageWriter
Expand Down Expand Up @@ -51,6 +50,7 @@ def __init__(self, name=None, tags=None, **kwargs):
self.mitiff_config = {}
self.translate_channel_name = {}
self.channel_order = {}
self.palette = False

def save_image(self):
raise NotImplementedError("save_image mitiff is not implemented.")
Expand All @@ -61,6 +61,8 @@ def save_dataset(self, dataset, filename=None, fill_value=None,

def _delayed_create(create_opts, dataset):
try:
if 'palette' in kwargs:
self.palette = kwargs['palette']
if 'platform_name' not in kwargs:
kwargs['platform_name'] = dataset.attrs['platform_name']
if 'name' not in kwargs:
Expand Down Expand Up @@ -160,6 +162,11 @@ def _make_channel_list(self, datasets, **kwargs):
channels.append(
ds.attrs['prerequisites'][ch][0])
break
elif self.palette:
if 'palette_channel_name' in kwargs:
channels.append(kwargs['palette_channel_name'].upper())
else:
LOG.error("Is palette but can not find palette_channel_name to name the dataset")
else:
for ch, ds in enumerate(datasets):
channels.append(ch + 1)
Expand Down Expand Up @@ -244,8 +251,8 @@ def _add_proj4_string(self, datasets, first_dataset):
if '+a=6378137.0 +b=6356752.31414' in proj4_string:
proj4_string = proj4_string.replace("+a=6378137.0 +b=6356752.31414",
"+ellps=WGS84")
if '+units=m' in proj4_string:
proj4_string = proj4_string.replace("+units=m", "+units=km")
if '+units=m' in proj4_string:
proj4_string = proj4_string.replace("+units=m", "+units=km")

if not any(datum in proj4_string for datum in ['datum', 'towgs84']):
proj4_string += ' +towgs84=0,0,0'
Expand All @@ -267,7 +274,20 @@ def _add_proj4_string(self, datasets, first_dataset):
proj4_string += ' +y_0=%.6f' % (
(-datasets.attrs['area'].area_extent[1] +
datasets.attrs['area'].pixel_size_y) + y_0)

elif '+x_0=0' in proj4_string and '+y_0=0' in proj4_string and isinstance(datasets, list):
proj4_string = proj4_string.replace("+x_0=0", '+x_0=%.6f' % (
(-first_dataset.attrs['area'].area_extent[0] +
first_dataset.attrs['area'].pixel_size_x) + x_0))
proj4_string = proj4_string.replace("+y_0=0", '+y_0=%.6f' % (
(-first_dataset.attrs['area'].area_extent[1] +
first_dataset.attrs['area'].pixel_size_y) + y_0))
elif '+x_0=0' in proj4_string and '+y_0=0' in proj4_string:
proj4_string = proj4_string.replace("+x_0=0", '+x_0=%.6f' % (
(-datasets.attrs['area'].area_extent[0] +
datasets.attrs['area'].pixel_size_x) + x_0))
proj4_string = proj4_string.replace("+y_0=0", '+y_0=%.6f' % (
(-datasets.attrs['area'].area_extent[1] +
datasets.attrs['area'].pixel_size_y) + y_0))
LOG.debug("proj4_string: %s", proj4_string)
proj4_string += '\n'

Expand Down Expand Up @@ -362,6 +382,17 @@ def _add_calibration_datasets(self, ch, datasets, reverse_offset, reverse_scale,
# http://stackoverflow.com/questions/1598579/rounding-decimals-with-new-python-format-function
return skip_calibration, _table_calibration, _reverse_offset, _reverse_scale, _decimals

def _add_palette_info(self, datasets, palette_unit, palette_description, **kwargs):
# mitiff key word for palette interpretion
_palette = '\n COLOR INFO:\n'
# mitiff info for the unit of the interpretion
_palette += ' {}\n'.format(palette_unit)
# The length of the palette description as needed by mitiff in DIANA
_palette += ' {}\n'.format(len(palette_description))
for desc in palette_description:
_palette += ' {}\n'.format(desc)
return _palette

def _add_calibration(self, channels, cns, datasets, **kwargs):

_table_calibration = ""
Expand Down Expand Up @@ -547,7 +578,11 @@ def _make_image_description(self, datasets, **kwargs):
else:
LOG.debug("Area extent: %s", datasets.attrs['area'].area_extent)

_image_description += self._add_calibration(channels, cns, datasets, **kwargs)
if self.palette:
LOG.debug("Doing palette image")
_image_description += self._add_palette_info(datasets, **kwargs)
else:
_image_description += self._add_calibration(channels, cns, datasets, **kwargs)

return _image_description

Expand All @@ -565,14 +600,56 @@ def _calibrate_data(self, dataset, calibration, min_val, max_val):
(float(max_val) - float(min_val))) * 255.
return _data.clip(0, 255)

def _save_as_palette(self, tif, datasets, **kwargs):
# MITIFF palette has only one data channel
if len(datasets.dims) == 2:
LOG.debug("Palette ok with only 2 dimensions. ie only x and y")
# 3 = Palette color. In this model, a color is described with a single component.
# The value of the component is used as an index into the red, green and blue curves
# in the ColorMap field to retrieve an RGB triplet that defines the color. When
# PhotometricInterpretation=3 is used, ColorMap must be present and SamplesPerPixel must be 1.
tif.SetField('PHOTOMETRIC', 3)

# As write_image can not save tiff image as palette, this has to be done basicly
# ie. all needed tags needs to be set.
tif.SetField('IMAGEWIDTH', datasets.sizes['x'])
tif.SetField('IMAGELENGTH', datasets.sizes['y'])
tif.SetField('BITSPERSAMPLE', 8)
tif.SetField('COMPRESSION', tif.get_tag_define('deflate'))

if 'palette_color_map' in kwargs:
tif.SetField('COLORMAP', kwargs['palette_color_map'])
else:
LOG.ERROR("In a mitiff palette image a color map must be provided: palette_color_map is missing.")

data_type = np.uint8
# Looks like we need to pass the data to writeencodedstrip as ctypes
tif.WriteEncodedStrip(0, np.ascontiguousarray(datasets.data.astype(data_type), data_type).ctypes.data,
datasets.sizes['x'] * datasets.sizes['y'])
tif.WriteDirectory()

def _save_as_enhanced(self, tif, datasets, **kwargs):
"""Save datasets as an enhanced RGB image"""
img = get_enhanced_image(datasets.squeeze(), enhance=self.enhancer)
if 'bands' in img.data.sizes and 'bands' not in datasets.sizes:
LOG.debug("Datasets without 'bands' become image with 'bands' due to enhancement.")
LOG.debug("Needs to regenerate mitiff image description")
image_description = self._make_image_description(img.data, **kwargs)
tif.SetField(IMAGEDESCRIPTION, (image_description).encode('utf-8'))
for i, band in enumerate(img.data['bands']):
chn = img.data.sel(bands=band)
data = chn.values.clip(0, 1) * 254. + 1
data = data.clip(0, 255)
tif.write_image(data.astype(np.uint8), compression='deflate')

def _save_datasets_as_mitiff(self, datasets, image_description,
gen_filename, **kwargs):
"""Put all togehter and save as a tiff file with the special tag
making it a mitiff file.
"""
from libtiff import TIFF

tif = TIFF.open(gen_filename, mode='w')
tif = TIFF.open(gen_filename, mode='wb')

tif.SetField(IMAGEDESCRIPTION, (image_description).encode('utf-8'))

Expand Down Expand Up @@ -618,18 +695,10 @@ def _save_datasets_as_mitiff(self, datasets, image_description,

tif.write_image(data.astype(np.uint8), compression='deflate')
break
elif self.palette:
LOG.debug("Saving dataset as palette.")
self._save_as_palette(tif, datasets, **kwargs)
else:
LOG.debug("Saving datasets as enhanced image")
img = get_enhanced_image(datasets.squeeze(), enhance=self.enhancer)
if 'bands' in img.data.sizes and 'bands' not in datasets.sizes:
LOG.debug("Datasets without 'bands' become image with 'bands' due to enhancement.")
LOG.debug("Needs to regenerate mitiff image description")
image_description = self._make_image_description(img.data, **kwargs)
tif.SetField(IMAGEDESCRIPTION, (image_description).encode('utf-8'))
for i, band in enumerate(img.data['bands']):
chn = img.data.sel(bands=band)
data = chn.values.clip(0, 1) * 254. + 1
data = data.clip(0, 255)
tif.write_image(data.astype(np.uint8), compression='deflate')

self._save_as_enhanced(tif, datasets, **kwargs)
del tif