diff --git a/satpy/etc/areas.def b/satpy/etc/areas.def index 7c6e0bbe82..d990539e31 100644 --- a/satpy/etc/areas.def +++ b/satpy/etc/areas.def @@ -44,14 +44,14 @@ REGION: iber { AREA_EXTENT: (-342379.698, 4432580.06, 723701.52, 5029648.75) }; -REGION: SouthAmerica { - NAME: South America - PCS_ID: geos0 - PCS_DEF: proj=geos,lon_0=0.0,a=6378169.00,b=6356583.80,h=35785831.0 - XSIZE: 3000 - YSIZE: 1200 - AREA_EXTENT: (-5570248.4773392612, -4263473.5610361192, -384719.90821206354, 1339786.2707295895) -}; +REGION: SouthAmerica { + NAME: South America + PCS_ID: geos0 + PCS_DEF: proj=geos,lon_0=0.0,a=6378169.00,b=6356583.80,h=35785831.0 + XSIZE: 3000 + YSIZE: 1200 + AREA_EXTENT: (-5570248.4773392612, -4263473.5610361192, -384719.90821206354, 1339786.2707295895) +}; REGION: brazil2 { NAME: brazil, platecarree @@ -72,12 +72,12 @@ REGION: sudeste { }; REGION: SouthAmerica_flat{ - NAME: South America flat + NAME: South America flat PCS_ID: eqc0 - PCS_DEF: proj='eqc', a=6378137.0, b=6378137.0, units='m' + PCS_DEF: proj='eqc', a=6378137.0, b=6378137.0, units='m' XSIZE: 1442 YSIZE: 1213 - AREA_EXTENT: (-8326322.8279089704,-4609377.085697311,-556597.45396636787, 1535833.8895192828 ) + AREA_EXTENT: (-8326322.8279089704,-4609377.085697311,-556597.45396636787, 1535833.8895192828 ) }; @@ -269,6 +269,24 @@ REGION: met09globeFull { AREA_EXTENT: (-5570248.4773392612, -5567248.074173444, 5567248.074173444, 5570248.4773392612) }; +REGION: seviri_0deg { + NAME: Full globe MSG image 0 degrees + PCS_ID: geos0 + PCS_DEF: proj=geos, lon_0=0, a=6378169.00, b=6356583.80, h=35785831.0 + XSIZE: 3712 + YSIZE: 3712 + AREA_EXTENT: (-5570248.4773392612, -5567248.074173444, 5567248.074173444, 5570248.4773392612) +}; + +REGION: seviri_iodc { + NAME: Full globe MSG image 41.5 degrees + PCS_ID: geos0 + PCS_DEF: proj=geos, lon_0=41.5, a=6378169.00, b=6356583.80, h=35785831.0 + XSIZE: 3712 + YSIZE: 3712 + AREA_EXTENT: (-5570248.4773392612, -5567248.074173444, 5567248.074173444, 5570248.4773392612) +}; + REGION: msg_resample_area { NAME: Full globe MSG image 20.75 degrees PCS_ID: geos0 @@ -567,7 +585,7 @@ REGION: eport1 { }; REGION: eport10 { - NAME: eport reduced resolution + NAME: eport reduced resolution PCS_ID: ps90n PCS_DEF: proj=stere,lat_0=90,lon_0=0,ellps=WGS84,units=m XSIZE: 1057 @@ -576,7 +594,7 @@ REGION: eport10 { }; REGION: eport4 { - NAME: eport reduced resolution + NAME: eport reduced resolution PCS_ID: ps90n PCS_DEF: proj=stere,lat_0=90,lon_0=0,ellps=WGS84,units=m XSIZE: 2642 @@ -585,7 +603,7 @@ REGION: eport4 { }; REGION: eport2 { - NAME: eport reduced resolution + NAME: eport reduced resolution PCS_ID: ps90n PCS_DEF: proj=stere,lat_0=90,lon_0=0,ellps=WGS84,units=m XSIZE: 5285 @@ -669,7 +687,7 @@ REGION: ease_nh { REGION: barents_sea { NAME: Barents and Greenland seas - PCS_ID: barents_sea + PCS_ID: barents_sea PCS_DEF: proj=stere,ellps=WGS84,lat_0=90,lon_0=40,lat_ts=75 XSIZE: 3000 YSIZE: 1700 @@ -809,7 +827,7 @@ REGION: northamerica_10km { # ------------------------------------------------------------------ -proj4_string = +proj4_string = REGION: romania { NAME: Romania - 3km @@ -894,4 +912,3 @@ REGION: robinson { YSIZE: 3296 AREA_EXTENT: (-2049911.5256036147, 5326895.725982913, 2049911.5256036168, 8625155.12857459) }; - diff --git a/satpy/etc/readers/native_msg.yaml b/satpy/etc/readers/native_msg.yaml index 8ae068ee7e..0945e566a0 100644 --- a/satpy/etc/readers/native_msg.yaml +++ b/satpy/etc/readers/native_msg.yaml @@ -8,9 +8,7 @@ reader: file_types: native_msg: file_reader: !!python/name:satpy.readers.native_msg.NativeMSGFileHandler '' - file_patterns: ['{satid:4s}-{instr:4s}-MSG{product_level:2d}-0100-NA-{processing_time1:%Y%m%d%H%M%S.%f}000Z-{processing_time2:%Y%m%d%H%M%S}-{order_id:d}.nat', - '{satid:4s}-{instr:4s}-MSG{product_level:2d}-0100-NA-{processing_time1:%Y%m%d%H%M%S.%f}000Z-NA.nat', - '{satid:4s}-{instr:4s}-MSG{product_level:2d}-0100-NA-{processing_time1:%Y%m%d%H%M%S.%f}000Z-{order_id:d}.nat'] + file_patterns: ['{satid:4s}-{instr:4s}-MSG{product_level:2d}-0100-NA-{processing_time1:%Y%m%d%H%M%S.%f}000Z-{order_id:s}.nat'] datasets: HRV: diff --git a/satpy/readers/native_msg.py b/satpy/readers/native_msg.py index 9bc8b61aff..5fa9ad84ee 100644 --- a/satpy/readers/native_msg.py +++ b/satpy/readers/native_msg.py @@ -7,7 +7,8 @@ # Adam Dybbroe # Ulrich Hamann -# Sauli Joro +# Sauli Joro +# Colin Duff # 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 @@ -63,7 +64,6 @@ def __init__(self, filename, filename_info, filetype_info): super(NativeMSGFileHandler, self).__init__(filename, filename_info, filetype_info) - self.filename = filename self.platform_name = None @@ -93,15 +93,13 @@ def _get_memmap(self): with open(self.filename) as fp_: - dt = self._get_filedtype() - - # Lazy reading: + dt = self._get_file_dtype() hdr_size = self.header.dtype.itemsize - return np.memmap(fp_, dtype=dt, shape=(self.data_len,), + return np.memmap(fp_, dtype=dt, shape=(self.mda['number_of_lines'],), offset=hdr_size, mode="r") - def _get_filedtype(self): + def _get_file_dtype(self): """Get the dtype of the file based on the actual available channels""" pkhrec = [ @@ -129,59 +127,56 @@ def get_lrec(cols): return lrec - visir_rec = get_lrec(self._cols_visir) - + visir_rec = get_lrec(int(self.mda['number_of_columns']*1.25)) number_of_lowres_channels = len( - [s for s in self._channel_index_list if not s == 'HRV']) + [s for s in self._channel_list if not s == 'HRV']) drec = [('visir', (visir_rec, number_of_lowres_channels))] if self.available_channels['HRV']: hrv_rec = get_lrec(int(self.mda['hrv_number_of_columns'] * 1.25)) drec.append(('hrv', (hrv_rec, 3))) - return np.dtype(drec) def _get_header(self): """Read the header info""" - hdrrec = Msg15NativeHeaderRecord().get() - hd_dt = np.dtype(hdrrec) - hd_dt = hd_dt.newbyteorder('>') - self.header = np.fromfile(self.filename, dtype=hd_dt, count=1) + hdr_rec = Msg15NativeHeaderRecord().get() + hdr_dtype = np.dtype(hdr_rec).newbyteorder('>') + self.header = np.fromfile(self.filename, dtype=hdr_dtype, count=1) + + data15hd = self.header['15_DATA_HEADER'] + sec15hd = self.header['15_SECONDARY_PRODUCT_HEADER'] # Set the list of available channels: self.available_channels = get_available_channels(self.header) - self._channel_index_list = [i for i in CHANNEL_NAMES.values() - if self.available_channels[i]] + self._channel_list = [i for i in CHANNEL_NAMES.values() + if self.available_channels[i]] - self.platform_id = self.header['15_DATA_HEADER'][ + self.platform_id = data15hd[ 'SatelliteStatus']['SatelliteDefinition']['SatelliteId'][0] self.platform_name = "Meteosat-" + SATNUM[self.platform_id] - ssp_lon = self.header['15_DATA_HEADER']['ImageDescription'][ - 'ProjectionDescription']['LongitudeOfSSP'][0] - self.mda = {} - equator_radius = self.header['15_DATA_HEADER']['GeometricProcessing'][ + + equator_radius = data15hd['GeometricProcessing'][ 'EarthModel']['EquatorialRadius'][0] * 1000. - north_polar_radius = self.header['15_DATA_HEADER'][ + north_polar_radius = data15hd[ 'GeometricProcessing']['EarthModel']['NorthPolarRadius'][0] * 1000. - south_polar_radius = self.header['15_DATA_HEADER'][ + south_polar_radius = data15hd[ 'GeometricProcessing']['EarthModel']['SouthPolarRadius'][0] * 1000. polar_radius = (north_polar_radius + south_polar_radius) * 0.5 + ssp_lon = data15hd['ImageDescription'][ + 'ProjectionDescription']['LongitudeOfSSP'][0] + self.mda['projection_parameters'] = {'a': equator_radius, 'b': polar_radius, 'h': 35785831.00, - 'SSP_longitude': ssp_lon} - - sec15hd = self.header['15_SECONDARY_PRODUCT_HEADER'] - numlines_visir = int(sec15hd['NumberLinesVISIR']['Value'][0]) - - self.mda['number_of_lines'] = numlines_visir + 'ssp_longitude': ssp_lon} west = int(sec15hd['WestColumnSelectedRectangle']['Value'][0]) east = int(sec15hd['EastColumnSelectedRectangle']['Value'][0]) north = int(sec15hd["NorthLineSelectedRectangle"]['Value'][0]) south = int(sec15hd["SouthLineSelectedRectangle"]['Value'][0]) + ncols_hrv_hdr = int(sec15hd['NumberColumnsHRV']['Value'][0]) # We suspect the UMARF will pad out any ROI colums that # arent divisible by 4 so here we work out how many pixels have @@ -191,45 +186,64 @@ def _get_header(self): if y < 4: # column has been padded with y pixels - self._cols_visir = int((west - east + 1 + y) * 1.25) + cols_visir = int((west - east + 1 + y) * 1.25) else: # no padding has occurred - self._cols_visir = int((west - east + 1) * 1.25) + cols_visir = int((west - east + 1) * 1.25) - self.mda['number_of_columns'] = int(self._cols_visir / 1.25) + # HRV Channel - check if an ROI + if (west - east) < 3711: + cols_hrv = int(np.ceil(ncols_hrv_hdr * 10.0 / 8)) # 6960 + else: + cols_hrv = int(np.ceil(5568 * 10.0 / 8)) # 6960 + + # self.mda should represent the 16bit dimensions not 10bit + self.mda['number_of_lines'] = int(sec15hd['NumberLinesVISIR']['Value'][0]) + self.mda['number_of_columns'] = int(cols_visir / 1.25) self.mda['hrv_number_of_lines'] = int(sec15hd["NumberLinesHRV"]['Value'][0]) - # The number of HRV columns seem to be correct in the UMARF header: - self.mda['hrv_number_of_columns'] = int(sec15hd["NumberColumnsHRV"]['Value'][0]) + self.mda['hrv_number_of_columns'] = int(cols_hrv / 1.25) + + # Check the calculated row,column dimensions against the header information: + ncols = self.mda['number_of_columns'] + ncols_hdr = int(sec15hd['NumberLinesVISIR']['Value'][0]) + + if ncols != ncols_hdr: + logger.warning( + "Number of VISIR columns from header and derived from data are not consistent!") + logger.warning("Number of columns read from header = %d", ncols_hdr) + logger.warning("Number of columns calculated from data = %d", ncols) + + ncols_hrv = self.mda['hrv_number_of_columns'] + + if ncols_hrv != ncols_hrv_hdr: + logger.warning( + "Number of HRV columns from header and derived from data are not consistent!") + logger.warning("Number of columns read from header = %d", ncols_hrv_hdr) + logger.warning("Number of columns calculated from data = %d", ncols_hrv) - coldir_step = self.header['15_DATA_HEADER']['ImageDescription'][ + column_step = data15hd['ImageDescription'][ "ReferenceGridVIS_IR"]["ColumnDirGridStep"][0] * 1000.0 - lindir_step = self.header['15_DATA_HEADER']['ImageDescription'][ + line_step = data15hd['ImageDescription'][ "ReferenceGridVIS_IR"]["LineDirGridStep"][0] * 1000.0 - # Check the calculated row,column dimensions against the header information: - numcols_cal = self.mda['number_of_columns'] - numcols_hd = self.header['15_DATA_HEADER']['ImageDescription']['ReferenceGridVIS_IR']['NumberOfColumns'][0] - if numcols_cal != numcols_hd: - logger.warning("Number of (non HRV band) columns from header and derived from data are not consistent!") - logger.warning("Number of columns read from header = %d", numcols_hd) - logger.warning("Number of columns calculated from data = %d", numcols_cal) - - area_extent = ((1856 - west - 0.5) * coldir_step, - (1856 - north + 0.5) * lindir_step, - (1856 - east + 0.5) * coldir_step, - (1856 - south + 1.5) * lindir_step) + # get corner points for lower left and upper right + # columns and lines + ll_c = (1856 - west - 0.5) * column_step + ll_l = (south - 1856 + 1.5) * line_step + ur_c = (1856 - east + 0.5) * column_step + ur_l = (north - 1856 + 0.5) * line_step - self.area_extent = area_extent + area_extent = (ll_c, ll_l, ur_c, ur_l) - self.data_len = numlines_visir + self.area_extent = area_extent def get_area_def(self, dsid): a = self.mda['projection_parameters']['a'] b = self.mda['projection_parameters']['b'] h = self.mda['projection_parameters']['h'] - lon_0 = self.mda['projection_parameters']['SSP_longitude'] + lon_0 = self.mda['projection_parameters']['ssp_longitude'] proj_dict = {'a': float(a), 'b': float(b), @@ -255,25 +269,27 @@ def get_area_def(self, dsid): self.area_extent) self.area = area + return area def get_dataset(self, dataset_id, info, xslice=slice(None), yslice=slice(None)): - channel_name = dataset_id.name + channel = dataset_id.name + channel_list = self._channel_list - if channel_name not in self._channel_index_list: - raise KeyError('Channel % s not available in the file' % channel_name) - elif channel_name not in ['HRV']: + if channel not in channel_list: + raise KeyError('Channel % s not available in the file' % channel) + elif channel not in ['HRV']: shape = (self.mda['number_of_lines'], self.mda['number_of_columns']) - ch_idn = self._channel_index_list.index(channel_name) # Check if there is only 1 channel in the list as a change # is needed in the arrray assignment ie channl id is not present - if len(self._channel_index_list) == 1: + if len(channel_list) == 1: raw = self.dask_array['visir']['line_data'] else: - raw = self.dask_array['visir']['line_data'][:, ch_idn, :] + i = channel_list.index(channel) + raw = self.dask_array['visir']['line_data'][:, i, :] data = mb.dec10216(raw.flatten()) data = da.flipud(da.fliplr((data.reshape(shape)))) @@ -301,63 +317,67 @@ def get_dataset(self, dataset_id, info, idx = range(2, shape[0], 3) data[idx, :] = data0 - res = xr.DataArray(data, dims=['y', 'x']).where(data != 0).astype(np.float32) + xarr = xr.DataArray(data, dims=['y', 'x']).where(data != 0).astype(np.float32) - if res is not None: - out = res + if xarr is None: + dataset = None else: - return None - - out = self.calibrate(out, dataset_id) - out.attrs['units'] = info['units'] - out.attrs['wavelength'] = info['wavelength'] - out.attrs['standard_name'] = info['standard_name'] - out.attrs['platform_name'] = self.platform_name - out.attrs['sensor'] = 'seviri' + dataset = self.calibrate(xarr, dataset_id) + dataset = xarr + dataset.attrs['units'] = info['units'] + dataset.attrs['wavelength'] = info['wavelength'] + dataset.attrs['standard_name'] = info['standard_name'] + dataset.attrs['platform_name'] = self.platform_name + dataset.attrs['sensor'] = 'seviri' - return out + return dataset def calibrate(self, data, dataset_id): """Calibrate the data.""" tic = datetime.now() + data15hdr = self.header['15_DATA_HEADER'] calibration = dataset_id.calibration - channel_name = dataset_id.name - channel_index = self._channel_index_list.index(channel_name) + channel = dataset_id.name + + # even though all the channels may not be present in the file, + # the header does have calibration coefficients for all the channels + # hence, this channel index needs to refer to full channel list + i = list(CHANNEL_NAMES.values()).index(channel) if calibration == 'counts': return if calibration in ['radiance', 'reflectance', 'brightness_temperature']: + # you cant apply GSICS values to the VIS channels + visual_channels = ['HRV', 'VIS006', 'VIS008', 'IR_016'] - calMode = 'NOMINAL' # determine the required calibration coefficients to use # for the Level 1.5 Header - # NB gsics doesnt apply to VIS channels so ignore them - if (channel_index > 2): - calMode = os.environ.get('CAL_MODE', 'NOMINAL') + calMode = os.environ.get('CAL_MODE', 'NOMINAL') - if (calMode.upper() != 'GSICS'): - coeffs = self.header['15_DATA_HEADER'][ + # NB GSICS doesn't have calibration coeffs for VIS channels + if (calMode.upper() != 'GSICS' or channel in visual_channels): + coeffs = data15hdr[ 'RadiometricProcessing']['Level15ImageCalibration'] - gain = coeffs['CalSlope'][0][channel_index] - offset = coeffs['CalOffset'][0][channel_index] + gain = coeffs['CalSlope'][0][i] + offset = coeffs['CalOffset'][0][i] else: - coeffs = self.header['15_DATA_HEADER'][ + coeffs = data15hdr[ 'RadiometricProcessing']['MPEFCalFeedback'] - gain = coeffs['GSICSCalCoeff'][0][channel_index] - offset = coeffs['GSICSOffsetCount'][0][channel_index] + gain = coeffs['GSICSCalCoeff'][0][i] + offset = coeffs['GSICSOffsetCount'][0][i] offset = offset * gain res = self._convert_to_radiance(data, gain, offset) if calibration == 'reflectance': - solar_irradiance = CALIB[self.platform_id][channel_name]["F"] + solar_irradiance = CALIB[self.platform_id][channel]["F"] res = self._vis_calibrate(res, solar_irradiance) elif calibration == 'brightness_temperature': - cal_type = self.header['15_DATA_HEADER']['ImageDescription'][ - 'Level15ImageProduction']['PlannedChanProcessing'][0][channel_index] - res = self._ir_calibrate(res, channel_name, cal_type) + cal_type = data15hdr['ImageDescription'][ + 'Level15ImageProduction']['PlannedChanProcessing'][0][i] + res = self._ir_calibrate(res, channel, cal_type) logger.debug("Calibration time " + str(datetime.now() - tic)) return res @@ -366,17 +386,18 @@ def calibrate(self, data, dataset_id): def get_available_channels(header): """Get the available channels from the header information""" + # Python2 and Python3 handle strings differently try: chlist_str = header['15_SECONDARY_PRODUCT_HEADER'][ 'SelectedBandIDs'][0][-1].strip().decode() except AttributeError: - # Strings have no deocde method in py3 + # Strings have no decode method in py3 chlist_str = header['15_SECONDARY_PRODUCT_HEADER'][ 'SelectedBandIDs'][0][-1].strip() retv = {} - for idx, chmark in zip(range(12), chlist_str): - retv[CHANNEL_NAMES[idx + 1]] = (chmark == 'X') + for idx, char in zip(range(12), chlist_str): + retv[CHANNEL_NAMES[idx + 1]] = (char == 'X') return retv diff --git a/satpy/readers/native_msg_hdr.py b/satpy/readers/native_msg_hdr.py index 691b4d92a3..13273011fd 100644 --- a/satpy/readers/native_msg_hdr.py +++ b/satpy/readers/native_msg_hdr.py @@ -22,11 +22,9 @@ # along with this program. If not, see . """Definition of Header Records for the MSG Level 1.5 data (hrit or native) - -.. warning:: - - `impf_configuration` in `L15DataHeaderRecord` class needs to be fixed! - + `impf_configuration` in `L15DataHeaderRecord` + Padding of 4 bytes has been added, believed to be down to Structure + padding by the compiler. """ @@ -34,7 +32,6 @@ class GSDTRecords(object): - """MSG Ground Segment Data Type records. Reference Document: MSG Ground Segment Design Specification (GSDS) @@ -98,8 +95,11 @@ class GSDTRecords(object): class Msg15NativeHeaderRecord(object): + """ + SEVIRI Level 1.5 header for native-format + """ - def get(self, umarf=True): + def get(self): record = [ ('15_MAIN_PRODUCT_HEADER', L15MainProductHeaderRecord().get()), @@ -107,7 +107,7 @@ def get(self, umarf=True): L15SecondaryProductHeaderRecord().get()), ('GP_PK_HEADER', GSDTRecords.gp_pk_header), ('GP_PK_SH1', GSDTRecords.gp_pk_sh1), - ('15_DATA_HEADER', L15DataHeaderRecord().get(umarf=umarf)) + ('15_DATA_HEADER', L15DataHeaderRecord().get()) ] return record @@ -122,7 +122,6 @@ class L15PhData(object): class L15MainProductHeaderRecord(object): - """ Reference Document: MSG Level 1.5 Native Format File Definition @@ -171,7 +170,6 @@ def get(self): class L15SecondaryProductHeaderRecord(object): - """ Reference Document: MSG Level 1.5 Native Format File Definition @@ -206,13 +204,12 @@ def get(self): class L15DataHeaderRecord(GSDTRecords): - """ Reference Document: MSG Level 1.5 Image Data Format Description """ - def get(self, umarf): + def get(self): record = [ ('15HeaderVersion', np.uint8), @@ -223,8 +220,7 @@ def get(self, umarf): ('RadiometricProcessing', self.radiometric_processing), ('GeometricProcessing', self.geometric_processing)] - if umarf: - record.append(('IMPFConfiguration', self.impf_configuration)) + record.append(('IMPFConfiguration', self.impf_configuration)) return record @@ -679,6 +675,6 @@ def impf_configuration(self): ('OverallConfiguration', overall_configuration), ('SUDetails', (su_details, 50)), ('WarmStartParams', warm_start_params), - ('Dummy', (np.void, 4))] # FIXME! + ('Padding', (np.void, 4))] return record