Source code for satpy.readers.generic_image

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2017 Satpy developers
#
# This file is part of satpy.
#
# satpy 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.
#
# satpy 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
# satpy.  If not, see <http://www.gnu.org/licenses/>.
# Author(s):

#   Lorenzo Clementi <lorenzo.clementi@meteoswiss.ch>
#   Panu Lahtinen <panu.lahtinen@fmi.fi>

# 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 <http://www.gnu.org/licenses/>.

"""
Reader for generic image (e.g. gif, png, jpg, tif, geotiff, ...).

Returns a dataset without calibration.  Includes coordinates if
available in the file (eg. geotiff).
"""

import logging

import xarray as xr
import dask.array as da
import numpy as np

from satpy.readers.file_handlers import BaseFileHandler
from satpy import CHUNK_SIZE

BANDS = {1: ['L'],
         2: ['L', 'A'],
         3: ['R', 'G', 'B'],
         4: ['R', 'G', 'B', 'A']}

logger = logging.getLogger(__name__)


[docs]class GenericImageFileHandler(BaseFileHandler): def __init__(self, filename, filename_info, filetype_info): super(GenericImageFileHandler, self).__init__( filename, filename_info, filetype_info) self.finfo = filename_info try: self.finfo['end_time'] = self.finfo['start_time'] except KeyError: pass self.finfo['filename'] = self.filename self.file_content = {} self.area = None self.read()
[docs] def read(self): """Read the image""" data = xr.open_rasterio(self.finfo["filename"], chunks=(1, CHUNK_SIZE, CHUNK_SIZE)) attrs = data.attrs.copy() # Create area definition if hasattr(data, 'crs'): self.area = self.get_geotiff_area_def(data.crs) # Rename to Satpy convention data = data.rename({'band': 'bands'}) # Rename bands to [R, G, B, A], or a subset of those data['bands'] = BANDS[data.bands.size] # Mask data if alpha channel is present try: data = mask_image_data(data) except ValueError as err: logger.warning(err) data.attrs = attrs self.file_content['image'] = data
[docs] def get_area_def(self, dsid): if self.area is None: raise NotImplementedError("No CRS information available from image") return self.area
[docs] def get_geotiff_area_def(self, crs): """Get extents from GeoTIFF file""" return get_geotiff_area_def(self.finfo['filename'], crs)
@property def start_time(self): return self.finfo['start_time'] @property def end_time(self): return self.finfo['end_time']
[docs] def get_dataset(self, key, info): """Get a dataset from the file.""" logger.debug("Reading %s.", key) return self.file_content[key.name]
[docs]def get_geotiff_area_def(filename, crs): """Read area definition from a geotiff.""" from osgeo import gdal from pyresample.geometry import AreaDefinition from pyresample.utils import proj4_str_to_dict fid = gdal.Open(filename) geo_transform = fid.GetGeoTransform() pcs_id = fid.GetProjection().split('"')[1] min_x = geo_transform[0] max_y = geo_transform[3] x_size = fid.RasterXSize y_size = fid.RasterYSize max_x = min_x + geo_transform[1] * x_size min_y = max_y + geo_transform[5] * y_size area_extent = [min_x, min_y, max_x, max_y] area_def = AreaDefinition('geotiff_area', pcs_id, pcs_id, proj4_str_to_dict(crs), x_size, y_size, area_extent) return area_def
[docs]def mask_image_data(data): """Mask image data if alpha channel is present.""" if data.bands.size in (2, 4): if not np.issubdtype(data.dtype, np.integer): raise ValueError("Only integer datatypes can be used as a mask.") mask = data.data[-1, :, :] == np.iinfo(data.dtype).min data = data.astype(np.float64) masked_data = da.stack([da.where(mask, np.nan, data.data[i, :, :]) for i in range(data.shape[0])]) data.data = masked_data data = data.sel(bands=BANDS[data.bands.size - 1]) return data