1# Copyright (c) 2015,2016,2017,2019 MetPy Developers. 2# Distributed under the terms of the BSD 3-Clause License. 3# SPDX-License-Identifier: BSD-3-Clause 4"""Tools to process GINI-formatted products.""" 5 6import contextlib 7from datetime import datetime 8from enum import Enum 9from io import BytesIO 10from itertools import repeat # noqa: I202 11import logging 12import re 13 14import numpy as np 15from xarray import Variable 16from xarray.backends.common import AbstractDataStore 17from xarray.core.utils import FrozenDict 18 19from ._tools import Bits, IOBuffer, NamedStruct, open_as_needed, zlib_decompress_all_frames 20from ..package_tools import Exporter 21 22exporter = Exporter(globals()) 23log = logging.getLogger(__name__) 24 25 26def _make_datetime(s): 27 r"""Convert 7 bytes from a GINI file to a `datetime` instance.""" 28 year, month, day, hour, minute, second, cs = s 29 if year < 70: 30 year += 100 31 return datetime(1900 + year, month, day, hour, minute, second, 10000 * cs) 32 33 34def _scaled_int(s): 35 r"""Convert a 3 byte string to a signed integer value.""" 36 # Get leftmost bit (sign) as 1 (if 0) or -1 (if 1) 37 sign = 1 - ((s[0] & 0x80) >> 6) 38 39 # Combine remaining bits 40 int_val = (((s[0] & 0x7f) << 16) | (s[1] << 8) | s[2]) 41 log.debug('Source: %s Int: %x Sign: %d', ' '.join(hex(c) for c in s), int_val, sign) 42 43 # Return scaled and with proper sign 44 return (sign * int_val) / 10000. 45 46 47def _name_lookup(names): 48 r"""Create an io helper to convert an integer to a named value.""" 49 mapper = dict(zip(range(len(names)), names)) 50 51 def lookup(val): 52 return mapper.get(val, 'UnknownValue') 53 return lookup 54 55 56class GiniProjection(Enum): 57 r"""Represents projection values in GINI files.""" 58 59 mercator = 1 60 lambert_conformal = 3 61 polar_stereographic = 5 62 63 64@exporter.export 65class GiniFile(AbstractDataStore): 66 """A class that handles reading the GINI format satellite images from the NWS. 67 68 This class attempts to decode every byte that is in a given GINI file. 69 70 Notes 71 ----- 72 The internal data structures that things are decoded into are subject to change. 73 74 """ 75 76 missing = 255 77 wmo_finder = re.compile('(T\\w{3}\\d{2})[\\s\\w\\d]+\\w*(\\w{3})\r\r\n') 78 79 crafts = ['Unknown', 'Unknown', 'Miscellaneous', 'JERS', 'ERS/QuikSCAT', 'POES/NPOESS', 80 'Composite', 'DMSP', 'GMS', 'METEOSAT', 'GOES-7', 'GOES-8', 'GOES-9', 81 'GOES-10', 'GOES-11', 'GOES-12', 'GOES-13', 'GOES-14', 'GOES-15', 'GOES-16'] 82 83 sectors = ['NH Composite', 'East CONUS', 'West CONUS', 'Alaska Regional', 84 'Alaska National', 'Hawaii Regional', 'Hawaii National', 'Puerto Rico Regional', 85 'Puerto Rico National', 'Supernational', 'NH Composite', 'Central CONUS', 86 'East Floater', 'West Floater', 'Central Floater', 'Polar Floater'] 87 88 channels = ['Unknown', 'Visible', 'IR (3.9 micron)', 'WV (6.5/6.7 micron)', 89 'IR (11 micron)', 'IR (12 micron)', 'IR (13 micron)', 'IR (1.3 micron)', 90 'Reserved', 'Reserved', 'Reserved', 'Reserved', 'Reserved', 'LI (Imager)', 91 'PW (Imager)', 'Surface Skin Temp (Imager)', 'LI (Sounder)', 'PW (Sounder)', 92 'Surface Skin Temp (Sounder)', 'CAPE', 'Land-sea Temp', 'WINDEX', 93 'Dry Microburst Potential Index', 'Microburst Day Potential Index', 94 'Convective Inhibition', 'Volcano Imagery', 'Scatterometer', 'Cloud Top', 95 'Cloud Amount', 'Rainfall Rate', 'Surface Wind Speed', 'Surface Wetness', 96 'Ice Concentration', 'Ice Type', 'Ice Edge', 'Cloud Water Content', 97 'Surface Type', 'Snow Indicator', 'Snow/Water Content', 'Volcano Imagery', 98 'Reserved', 'Sounder (14.71 micron)', 'Sounder (14.37 micron)', 99 'Sounder (14.06 micron)', 'Sounder (13.64 micron)', 'Sounder (13.37 micron)', 100 'Sounder (12.66 micron)', 'Sounder (12.02 micron)', 'Sounder (11.03 micron)', 101 'Sounder (9.71 micron)', 'Sounder (7.43 micron)', 'Sounder (7.02 micron)', 102 'Sounder (6.51 micron)', 'Sounder (4.57 micron)', 'Sounder (4.52 micron)', 103 'Sounder (4.45 micron)', 'Sounder (4.13 micron)', 'Sounder (3.98 micron)', 104 # Percent Normal TPW found empirically from Service Change Notice 20-03 105 'Sounder (3.74 micron)', 'Sounder (Visible)', 'Percent Normal TPW'] 106 107 prod_desc_fmt = NamedStruct([('source', 'b'), 108 ('creating_entity', 'b', _name_lookup(crafts)), 109 ('sector_id', 'b', _name_lookup(sectors)), 110 ('channel', 'b', _name_lookup(channels)), 111 ('num_records', 'H'), ('record_len', 'H'), 112 ('datetime', '7s', _make_datetime), 113 ('projection', 'b', GiniProjection), ('nx', 'H'), ('ny', 'H'), 114 ('la1', '3s', _scaled_int), ('lo1', '3s', _scaled_int) 115 ], '>', 'ProdDescStart') 116 117 lc_ps_fmt = NamedStruct([('reserved', 'b'), ('lov', '3s', _scaled_int), 118 ('dx', '3s', _scaled_int), ('dy', '3s', _scaled_int), 119 ('proj_center', 'b')], '>', 'LambertOrPolarProjection') 120 121 mercator_fmt = NamedStruct([('resolution', 'b'), ('la2', '3s', _scaled_int), 122 ('lo2', '3s', _scaled_int), ('di', 'H'), ('dj', 'H') 123 ], '>', 'MercatorProjection') 124 125 prod_desc2_fmt = NamedStruct([('scanning_mode', 'b', Bits(3)), 126 ('lat_in', '3s', _scaled_int), ('resolution', 'b'), 127 ('compression', 'b'), ('version', 'b'), ('pdb_size', 'H'), 128 ('nav_cal', 'b')], '>', 'ProdDescEnd') 129 130 nav_fmt = NamedStruct([('sat_lat', '3s', _scaled_int), ('sat_lon', '3s', _scaled_int), 131 ('sat_height', 'H'), ('ur_lat', '3s', _scaled_int), 132 ('ur_lon', '3s', _scaled_int)], '>', 'Navigation') 133 134 def __init__(self, filename): 135 r"""Create an instance of `GiniFile`. 136 137 Parameters 138 ---------- 139 filename : str or file-like object 140 If str, the name of the file to be opened. Gzip-ed files are 141 recognized with the extension ``'.gz'``, as are bzip2-ed files with 142 the extension ``'.bz2'`` If `filename` is a file-like object, 143 this will be read from directly. 144 145 """ 146 fobj = open_as_needed(filename) 147 148 # Just read in the entire set of data at once 149 with contextlib.closing(fobj): 150 self._buffer = IOBuffer.fromfile(fobj) 151 152 # Pop off the WMO header if we find it 153 self.wmo_code = '' 154 self._process_wmo_header() 155 log.debug('First wmo code: %s', self.wmo_code) 156 157 # Decompress the data if necessary, and if so, pop off new header 158 log.debug('Length before decompression: %s', len(self._buffer)) 159 self._buffer = IOBuffer(self._buffer.read_func(zlib_decompress_all_frames)) 160 log.debug('Length after decompression: %s', len(self._buffer)) 161 162 # Process WMO header inside compressed data if necessary 163 self._process_wmo_header() 164 log.debug('2nd wmo code: %s', self.wmo_code) 165 166 # Read product description start 167 start = self._buffer.set_mark() 168 169 #: :desc: Decoded first section of product description block 170 #: :type: namedtuple 171 self.prod_desc = self._buffer.read_struct(self.prod_desc_fmt) 172 log.debug(self.prod_desc) 173 174 #: :desc: Decoded geographic projection information 175 #: :type: namedtuple 176 self.proj_info = None 177 178 # Handle projection-dependent parts 179 if self.prod_desc.projection in (GiniProjection.lambert_conformal, 180 GiniProjection.polar_stereographic): 181 self.proj_info = self._buffer.read_struct(self.lc_ps_fmt) 182 elif self.prod_desc.projection == GiniProjection.mercator: 183 self.proj_info = self._buffer.read_struct(self.mercator_fmt) 184 else: 185 log.warning('Unknown projection: %d', self.prod_desc.projection) 186 log.debug(self.proj_info) 187 188 # Read the rest of the guaranteed product description block (PDB) 189 #: :desc: Decoded second section of product description block 190 #: :type: namedtuple 191 self.prod_desc2 = self._buffer.read_struct(self.prod_desc2_fmt) 192 log.debug(self.prod_desc2) 193 194 if self.prod_desc2.nav_cal not in (0, -128): # TODO: See how GEMPAK/MCIDAS parses 195 # Only warn if there actually seems to be useful navigation data 196 if self._buffer.get_next(self.nav_fmt.size) != b'\x00' * self.nav_fmt.size: 197 log.warning('Navigation/Calibration unhandled: %d', self.prod_desc2.nav_cal) 198 if self.prod_desc2.nav_cal in (1, 2): 199 self.navigation = self._buffer.read_struct(self.nav_fmt) 200 log.debug(self.navigation) 201 202 # Catch bad PDB with size set to 0 203 if self.prod_desc2.pdb_size == 0: 204 log.warning('Adjusting bad PDB size from 0 to 512.') 205 self.prod_desc2 = self.prod_desc2._replace(pdb_size=512) 206 207 # Jump past the remaining empty bytes in the product description block 208 self._buffer.jump_to(start, self.prod_desc2.pdb_size) 209 210 # Read the actual raster--unless it's PNG compressed, in which case that happens later 211 blob = self._buffer.read(self.prod_desc.num_records * self.prod_desc.record_len) 212 213 # Check for end marker 214 end = self._buffer.read(self.prod_desc.record_len) 215 if end != b''.join(repeat(b'\xff\x00', self.prod_desc.record_len // 2)): 216 log.warning('End marker not as expected: %s', end) 217 218 # Check to ensure that we processed all of the data 219 if not self._buffer.at_end(): 220 if not blob: 221 log.debug('No data read yet, trying to decompress remaining data as an image.') 222 from matplotlib.image import imread 223 blob = (imread(BytesIO(self._buffer.read())) * 255).astype('uint8') 224 else: 225 log.warning('Leftover unprocessed data beyond EOF marker: %s', 226 self._buffer.get_next(10)) 227 228 self.data = np.array(blob).reshape((self.prod_desc.ny, 229 self.prod_desc.nx)) 230 231 def _process_wmo_header(self): 232 """Read off the WMO header from the file, if necessary.""" 233 data = self._buffer.get_next(64).decode('utf-8', 'ignore') 234 match = self.wmo_finder.search(data) 235 if match: 236 self.wmo_code = match.groups()[0] 237 self.siteID = match.groups()[-1] 238 self._buffer.skip(match.end()) 239 240 def __str__(self): 241 """Return a string representation of the product.""" 242 parts = [self.__class__.__name__ + ': {0.creating_entity} {0.sector_id} {0.channel}', 243 'Time: {0.datetime}', 'Size: {0.ny}x{0.nx}', 244 'Projection: {0.projection.name}', 245 'Lower Left Corner (Lon, Lat): ({0.lo1}, {0.la1})', 246 'Resolution: {1.resolution}km'] 247 return '\n\t'.join(parts).format(self.prod_desc, self.prod_desc2) 248 249 def _make_proj_var(self): 250 proj_info = self.proj_info 251 prod_desc2 = self.prod_desc2 252 attrs = {'earth_radius': 6371200.0} 253 if self.prod_desc.projection == GiniProjection.lambert_conformal: 254 attrs['grid_mapping_name'] = 'lambert_conformal_conic' 255 attrs['standard_parallel'] = prod_desc2.lat_in 256 attrs['longitude_of_central_meridian'] = proj_info.lov 257 attrs['latitude_of_projection_origin'] = prod_desc2.lat_in 258 elif self.prod_desc.projection == GiniProjection.polar_stereographic: 259 attrs['grid_mapping_name'] = 'polar_stereographic' 260 attrs['straight_vertical_longitude_from_pole'] = proj_info.lov 261 attrs['latitude_of_projection_origin'] = -90 if proj_info.proj_center else 90 262 attrs['standard_parallel'] = 60.0 # See Note 2 for Table 4.4A in ICD 263 elif self.prod_desc.projection == GiniProjection.mercator: 264 attrs['grid_mapping_name'] = 'mercator' 265 attrs['longitude_of_projection_origin'] = self.prod_desc.lo1 266 attrs['latitude_of_projection_origin'] = self.prod_desc.la1 267 attrs['standard_parallel'] = prod_desc2.lat_in 268 else: 269 raise NotImplementedError( 270 f'Unhandled GINI Projection: {self.prod_desc.projection}') 271 272 return 'projection', Variable((), 0, attrs) 273 274 def _make_time_var(self): 275 base_time = self.prod_desc.datetime.replace(hour=0, minute=0, second=0, microsecond=0) 276 offset = self.prod_desc.datetime - base_time 277 time_var = Variable((), data=offset.seconds + offset.microseconds / 1e6, 278 attrs={'units': 'seconds since ' + base_time.isoformat()}) 279 280 return 'time', time_var 281 282 def _get_proj_and_res(self): 283 import pyproj 284 285 proj_info = self.proj_info 286 prod_desc2 = self.prod_desc2 287 288 kwargs = {'a': 6371200.0, 'b': 6371200.0} 289 if self.prod_desc.projection == GiniProjection.lambert_conformal: 290 kwargs['proj'] = 'lcc' 291 kwargs['lat_0'] = prod_desc2.lat_in 292 kwargs['lon_0'] = proj_info.lov 293 kwargs['lat_1'] = prod_desc2.lat_in 294 kwargs['lat_2'] = prod_desc2.lat_in 295 dx, dy = proj_info.dx, proj_info.dy 296 elif self.prod_desc.projection == GiniProjection.polar_stereographic: 297 kwargs['proj'] = 'stere' 298 kwargs['lon_0'] = proj_info.lov 299 kwargs['lat_0'] = -90 if proj_info.proj_center else 90 300 kwargs['lat_ts'] = 60.0 # See Note 2 for Table 4.4A in ICD 301 kwargs['x_0'] = False # Easting 302 kwargs['y_0'] = False # Northing 303 dx, dy = proj_info.dx, proj_info.dy 304 elif self.prod_desc.projection == GiniProjection.mercator: 305 kwargs['proj'] = 'merc' 306 kwargs['lat_0'] = self.prod_desc.la1 307 kwargs['lon_0'] = self.prod_desc.lo1 308 kwargs['lat_ts'] = prod_desc2.lat_in 309 kwargs['x_0'] = False # Easting 310 kwargs['y_0'] = False # Northing 311 dx, dy = prod_desc2.resolution, prod_desc2.resolution 312 313 return pyproj.Proj(**kwargs), dx, dy 314 315 def _make_coord_vars(self): 316 proj, dx, dy = self._get_proj_and_res() 317 318 # Get projected location of lower left point 319 x0, y0 = proj(self.prod_desc.lo1, self.prod_desc.la1) 320 321 # Coordinate variable for x 322 xlocs = x0 + np.arange(self.prod_desc.nx) * (1000. * dx) 323 attrs = {'units': 'm', 'long_name': 'x coordinate of projection', 324 'standard_name': 'projection_x_coordinate'} 325 x_var = Variable(('x',), xlocs, attrs) 326 327 # Now y--Need to flip y because we calculated from the lower left corner, 328 # but the raster data is stored with top row first. 329 ylocs = (y0 + np.arange(self.prod_desc.ny) * (1000. * dy))[::-1] 330 attrs = {'units': 'm', 'long_name': 'y coordinate of projection', 331 'standard_name': 'projection_y_coordinate'} 332 y_var = Variable(('y',), ylocs, attrs) 333 334 # Get the two-D lon,lat grid as well 335 x, y = np.meshgrid(xlocs, ylocs) 336 lon, lat = proj(x, y, inverse=True) 337 338 lon_var = Variable(('y', 'x'), data=lon, 339 attrs={'long_name': 'longitude', 'units': 'degrees_east'}) 340 lat_var = Variable(('y', 'x'), data=lat, 341 attrs={'long_name': 'latitude', 'units': 'degrees_north'}) 342 343 return [('x', x_var), ('y', y_var), ('lon', lon_var), ('lat', lat_var)] 344 345 def get_variables(self): 346 """Get all variables in the file. 347 348 This is used by `xarray.open_dataset`. 349 350 """ 351 variables = [self._make_time_var()] 352 proj_var_name, proj_var = self._make_proj_var() 353 variables.append((proj_var_name, proj_var)) 354 variables.extend(self._make_coord_vars()) 355 356 # Now the data 357 name = self.prod_desc.channel 358 if '(' in name: 359 name = name.split('(')[0].rstrip() 360 361 missing_val = self.missing 362 attrs = {'long_name': self.prod_desc.channel, 'missing_value': missing_val, 363 'coordinates': 'y x time', 'grid_mapping': proj_var_name} 364 data_var = Variable(('y', 'x'), 365 data=np.ma.array(self.data, 366 mask=self.data == missing_val), 367 attrs=attrs) 368 variables.append((name, data_var)) 369 370 return FrozenDict(variables) 371 372 def get_attrs(self): 373 """Get the global attributes. 374 375 This is used by `xarray.open_dataset`. 376 377 """ 378 return FrozenDict(satellite=self.prod_desc.creating_entity, 379 sector=self.prod_desc.sector_id) 380