Source code for satpy.readers.li_l2

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2016 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/>.
"""Interface to MTG-LI L2 product NetCDF files

The reader is based on preliminary test data provided by EUMETSAT.
The data description is described in the
"LI L2 Product User Guide [LIL2PUG] Draft version" documentation.

"""
import h5netcdf
import logging
import numpy as np
from datetime import datetime
from pyresample import geometry
from satpy.readers.file_handlers import BaseFileHandler
# FIXME: This is not xarray/dask compatible
# TODO: Once migrated to xarray/dask, remove ignored path in setup.cfg
from satpy.dataset import Dataset

logger = logging.getLogger(__name__)


[docs]class LIFileHandler(BaseFileHandler): """MTG LI File Reader.""" def __init__(self, filename, filename_info, filetype_info): super(LIFileHandler, self).__init__(filename, filename_info, filetype_info) self.nc = h5netcdf.File(self.filename, 'r') # Get grid dimensions from file refdim = self.nc['grid_position'][:] # Get number of lines and columns self.nlines = int(refdim[2]) self.ncols = int(refdim[3]) self.cache = {} logger.debug('Dimension : {}'.format(refdim)) logger.debug('Row/Cols: {} / {}'.format(self.nlines, self.ncols)) logger.debug('Reading: {}'.format(self.filename)) logger.debug('Start: {}'.format(self.start_time)) logger.debug('End: {}'.format(self.end_time)) @property def start_time(self): return datetime.strptime(self.nc.attrs['sensing_start'], '%Y%m%d%H%M%S') @property def end_time(self): return datetime.strptime(self.nc.attrs['end_time'], '%Y%m%d%H%M%S')
[docs] def get_dataset(self, key, info=None, out=None): """Load a dataset """ if key in self.cache: return self.cache[key] # Type dictionary typedict = {"af": "flash_accumulation", "afa": "accumulated_flash_area", "afr": "flash_radiance", "lgr": "radiance", "lef": "radiance", "lfl": "radiance"} # Get lightning data out of NetCDF container logger.debug("Key: {}".format(key.name)) # Create reference grid grid = np.full((self.nlines, self.ncols), np.NaN) # Get product values values = self.nc[typedict[key.name]] rows = self.nc['row'] cols = self.nc['column'] logger.debug('[ Number of values ] : {}'.format((len(values)))) logger.debug('[Min/Max] : <{}> / <{}>'.format(np.min(values), np.max(values))) # Convert xy coordinates to flatten indices ids = np.ravel_multi_index([rows, cols], grid.shape) # Replace NaN values with data np.put(grid, ids, values) # Correct for bottom left origin in LI row/column indices. rotgrid = np.flipud(grid) # Rotate the grid by 90 degree clockwise rotgrid = np.rot90(rotgrid, 3) logger.warning("LI data has been rotated to fit to reference grid. \ Works only for test dataset") # Mask invalid values ds = np.ma.masked_where(np.isnan(rotgrid), rotgrid) # Create dataset object out.data[:] = np.ma.getdata(ds) out.mask[:] = np.ma.getmask(ds) out.info.update(key.to_dict()) return out
[docs] def get_area_def(self, key, info=None): """Create AreaDefinition for specified product. Projection information are hard coded for 0 degree geos projection Test dataset doesn't provide the values in the file container. Only fill values are inserted. """ # TODO Get projection information from input file a = 6378169. h = 35785831. b = 6356583.8 lon_0 = 0. # area_extent = (-5432229.9317116784, -5429229.5285458621, # 5429229.5285458621, 5432229.9317116784) area_extent = (-5570248.4773392612, -5567248.074173444, 5567248.074173444, 5570248.4773392612) proj_dict = {'a': float(a), 'b': float(b), 'lon_0': float(lon_0), 'h': float(h), 'proj': 'geos', 'units': 'm'} area = geometry.AreaDefinition( 'LI_area_name', "LI area", 'geosli', proj_dict, self.ncols, self.nlines, area_extent) self.area = area logger.debug("Dataset area definition: \n {}".format(area)) return area