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

Calculate gradient and offset values from satpy stretching data #15

Draft
wants to merge 6 commits into
base: master
Choose a base branch
from
Draft

Calculate gradient and offset values from satpy stretching data #15

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
ef5fa09
introduced reading scaling values range from satpy processing
vgiuffrida May 15, 2019
ce1ae14
use satpy stretch values for slope data
vgiuffrida Aug 10, 2019
40f0d6f
remove comment
vgiuffrida Aug 10, 2019
9b70664
fix use_value_range values reading
vgiuffrida Aug 14, 2019
c37aea7
remove spaces
vgiuffrida Aug 14, 2019
d300943
remove comment
vgiuffrida Aug 16, 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
124 changes: 115 additions & 9 deletions pyninjotiff/ninjotiff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -316,7 +316,7 @@ def _get_satellite_altitude(filename):


def _finalize(img, dtype=np.uint8, value_range_measurement_unit=None,
data_is_scaled_01=True, fill_value=None):
data_is_scaled_01=True, fill_value=None, use_value_range=True):
"""Finalize a mpop GeoImage for Ninjo. Specialy take care of phycical scale
and offset.

Expand Down Expand Up @@ -349,20 +349,78 @@ def _finalize(img, dtype=np.uint8, value_range_measurement_unit=None,
They still have to be tested thoroughly.
"""

if img.mode == 'L':
# use_value_range mode: calculate slope and offset from satpy stretching values
if use_value_range:

if fill_value is not None:
log.debug("Forcing fill value to %s", fill_value)

data = img.data

fill_value = fill_value if fill_value is not None else np.iinfo(dtype).min
log.debug("Entry xarray: %.2f, %.2f, %.2f" %
(data.min(), data.mean(), data.max() ))

if (value_range_measurement_unit is not None and
len(value_range_measurement_unit)>=2):
if len(value_range_measurement_unit)<4:
value_range_measurement_unit.append(value_range_measurement_unit[0])
value_range_measurement_unit.append(value_range_measurement_unit[1])
value_range_measurement_unit[0] = 0.
value_range_measurement_unit[1] = 1.
data = data.clip(value_range_measurement_unit[0],value_range_measurement_unit[1])
log.debug("After clipping: %.2f, %.2f, %.2f" %
(data.min(), data.mean(), data.max()))

scale_val = ((value_range_measurement_unit[1] - value_range_measurement_unit[0]) /
(np.iinfo(dtype).max))
scale_val = scale_val or 1
offset_val = value_range_measurement_unit[0]

scale = ((value_range_measurement_unit[3] - value_range_measurement_unit[2]) /
(np.iinfo(dtype).max))
scale = scale or 1
offset = value_range_measurement_unit[2]
else:
scale_val = scale = 1
scale_val = 1. / np.iinfo(dtype).max
offset_val = offset = data.min()

log.debug("Before scaling: %.2f, %.2f, %.2f" %
(data.min(), data.mean(), data.max()))


scale_fill = ((np.iinfo(dtype).max) / (np.iinfo(dtype).max + 1.0))
data *= scale_fill

data = (1 + ((data-offset_val) / scale_val)).astype(dtype)

log.debug("After scaling: %.2f, %.2f, %.2f" %
(data.min(), data.mean(), data.max()))

return data[0], scale, offset, fill_value


if img.mode == 'L' and not use_value_range:
# PFE: mpop.satout.cfscene
if isinstance(img, np.ma.MaskedArray):
data = img.channels[0]
else :
# TODO: check what is the corret fill value for NinJo!
if fill_value is not None:
log.debug("Forcing fill value to %s", fill_value)
non_band_dims = tuple(x for x in img.data.dims if x != 'bands')
test_chn_max = img.data[0].max()
test_chn_min = img.data[0].min()
test_scale = ( ((test_chn_max - test_chn_min) * 255) / (np.iinfo(dtype).max - 1.0) )
data = img.finalize(dtype=dtype, fill_value=fill_value)
test_chn_max2 = data[0].max()
test_chn_min2 = data[0].min()
# Go back to the masked_array for compatibility
# with the following part of the code.
data = data[0].to_masked_array()

fill_value = fill_value if fill_value is not None else np.iinfo(dtype).min
fill_value = fill_value if fill_value is not None else np.iinfo(dtype).min

log.debug("Before scaling: %.2f, %.2f, %.2f" %
(data.min(), data.mean(), data.max()))
Expand All @@ -385,8 +443,15 @@ def _finalize(img, dtype=np.uint8, value_range_measurement_unit=None,
img = deepcopy(img)
img.channels[0] *= scale_fill_value

#test
#data = img.data
#data[0] *= scale_fill_value

img.channels[0] += 1 / (np.iinfo(dtype).max + 1.0)

#test
#data[0] += 1 / (np.iinfo(dtype).max + 1.0)

channels, fill_value = img._finalize(dtype)
data = channels[0]

Expand All @@ -403,7 +468,6 @@ def _finalize(img, dtype=np.uint8, value_range_measurement_unit=None,

if fill_value is None:
fill_value = 0

else:
if value_range_measurement_unit:
data.clip(value_range_measurement_unit[0],
Expand All @@ -413,8 +477,13 @@ def _finalize(img, dtype=np.uint8, value_range_measurement_unit=None,
log.debug("Scaling, using value range %.2f - %.2f" %
(value_range_measurement_unit[0], value_range_measurement_unit[1]))
else:

chn_max = data.max()
chn_min = data.min()

#test
#chn_max = data[0].max()
#chn_min = data[0].min()
log.debug("Doing auto scaling")

# Make room for transparent pixel.
Expand Down Expand Up @@ -545,12 +614,48 @@ def save(img, filename, ninjo_product_name=None, writer_options=None,
if 'fill_value' in kwargs and kwargs['fill_value'] != None:
fill_value = int(kwargs['fill_value'])

value_range_measurement_unit = None

# If use_value_range is set, use stretch values from satpy
# to calculate gradien and offset for ninjotiff tags
use_range = False if 'use_value_range' not in kwargs else kwargs['use_value_range']

# If new mode for pyninjotiff
if 'output_info' in img.data.attrs and use_range:
min_value = None
max_value = None
if isinstance(img.data.attrs['output_info']['max_value'],float):
min_value = img.data.attrs['output_info']['min_value']
max_value = img.data.attrs['output_info']['max_value']
else:
try:
max_value = img.data.attrs['output_info']['max_value'][0]
min_value = img.data.attrs['output_info']['min_value'][0]
except IndexError:
max_value = img.data.attrs['output_info']['max_value'][()]
min_value = img.data.attrs['output_info']['min_value'][()]
except KeyError:
log.debug("output_info data not found")
if min_value is not None and max_value is not None:
# do we need to convert kelvin to celcius
value_range_measurement_unit = [min_value, max_value]
convert = 0.
if ('physic_value' in img.data.attrs['output_info'] and
'physic_unit' in img.data.attrs['output_info']):
# If physic_value and physic_unit are set, adjust min and max values using the correct unit
if img.data.attrs['output_info']['physic_value'] == "T" :
if img.data.attrs['output_info']['physic_unit'] == "KELVIN" and kwargs["physic_unit"] == "CELSIUS":
convert = -273
if img.data.attrs['output_info']['physic_unit'] == "CELSIUS" and kwargs["physic_unit"] == "KELVIN":
convert = 273
value_range_measurement_unit = [float(min_value)+convert,float(max_value)+convert]

try:
value_range_measurement_unit = (float(kwargs["ch_min_measurement_unit"]),
float(kwargs["ch_max_measurement_unit"]))
value_range_measurement_unit = [float(kwargs["ch_min_measurement_unit"]),
float(kwargs["ch_max_measurement_unit"])]
except KeyError:
value_range_measurement_unit = None
log.debug("not forced ch_min and max in ninjo config file")


data_is_scaled_01 = bool(kwargs.get("data_is_scaled_01", True))

Expand All @@ -561,7 +666,8 @@ def save(img, filename, ninjo_product_name=None, writer_options=None,
dtype=dtype,
data_is_scaled_01=data_is_scaled_01,
value_range_measurement_unit=value_range_measurement_unit,
fill_value=fill_value,)
fill_value=fill_value,
use_value_range=use_range)

if isinstance(img, np.ma.MaskedArray):
area_def = img.info['area']
Expand Down Expand Up @@ -686,7 +792,7 @@ def write(image_data, output_fn, area_def, product_name=None, **kwargs):
options['ref_lat2'] = 0
if 'lon_0' in area_def.proj_dict:
options['central_meridian'] = area_def.proj_dict['lon_0']


a,b = proj4_radius_parameters(area_def.proj_dict)
options['radius_a'] = a
Expand Down