Save Datasets using SCMI

[WxmanJ]

I am having an issue using the beta scmi writer. I spoke with David and he mentioned several of the variables are not needed; however, the writer seems to give warnings/errors unless I add the platform, sensor, source_name, and others. Upon entering all requested variables, I am still getting the KeyError: 'platform'.

Here is my script.

from satpy.scene import Scene
from datetime import datetime, timedelta
from netCDF4 import Dataset
from glob import glob
import numpy as np
from pyproj import Proj
from satpy.writers import Writer, DecisionTree, scmi
from pyresample.geometry import AreaDefinition

FILES = glob("/root/*.DAT")

global_scene = Scene(reader='ahi_hsd',filenames=FILES)
global_scene.load([0.47])
print(global_scene[0.47])

global_scene.save_datasets(writer='scmi', sector_id="Pacific", source_name="AWIPS", platform_name="Himawari", filename="AWIPS_AII_HIMAWARI_B01.nc")

Here is the traceback I receive using debug.

[INFO: 2018-03-14 23:57:08 : satpy.writers.scmi] Writing product B01 to AWIPS SCMI NetCDF file
[ERROR: 2018-03-14 23:57:08 : satpy.writers.scmi] Could not get information on dataset from backend configuration file
Traceback (most recent call last):
  File "test3.py", line 23, in <module>
    global_scene.save_datasets(writer='scmi', sector_id="Pacific", sensor="AHI", source_name="AWIPS", platform_name="Himawari", filename="AWIPS_AII_HIMAWARI_B01.nc")
  File "/usr/lib/python2.7/site-packages/satpy/scene.py", line 672, in save_datasets
    writer.save_datasets(datasets, **kwargs)
  File "/usr/lib/python2.7/site-packages/satpy/writers/scmi.py", line 848, in save_datasets
    pkwargs['awips_info'] = self._get_awips_info(ds_info, source_name=source_name)
  File "/usr/lib/python2.7/site-packages/satpy/writers/scmi.py", line 806, in _get_awips_info
    def_ce = "{}-{}".format(ds_info["platform"].upper(), ds_info["sensor"].upper())
KeyError: 'platform'

[WxmanJ]

I was able to resolve the issue by defining the variable in scmi.py. I opted to call it "platform_name" for my test (below). Thereafter, the scmi writer ran smoothly and created the image. I copied the image to AWIPS and it does indeed work as expected.

def save_datasets(self, datasets, sector_id=None,
                  source_name=None, filename=None,
                  tile_count=(1, 1), tile_size=None,
                  lettered_grid=False, num_subtiles=None,
	          platform_name=None,
                  **kwargs):

[djhoese]

@WxmanJ Could you please undo your changes to scmi.py and we'll test a few things out. Please remove the sensor="AHI" from your call to save_datasets as well.

The way the SCMI writer (and all of satpy works) is by using the metadata attached to each dataset to figure out what that dataset is and what to do with it. In the version of satpy you are working with all of a dataset's metadata is kept in a ds.info dictionary. The ds_info in your error message above is this ds.info dictionary. It is complaining because there doesn't seem to be a 'platform' key in your ds.info. Could you please tell me what print(global_scene[0.47].info) shows?

[WxmanJ]

Adding "sensor" was actually a mistyping on my part. The script I am using does not actually contain that -- it was in the original script, but not the modified one I am using now. I edited the issue and removed it. I also removed my changes to scmi.py.

Adding the print statement requested, this is what is output.

{'polarization': None, 'sensor': 'ahi', 'start_time': datetime.datetime(2018, 3, 15, 6, 0, 21, 136857), 'name': 'B01', 'satellite_longitude': 140.64899219645, 'area': Area ID: some_area_name
Description: On-the-fly area
Projection ID: geosh8
Projection: {'a': '6378137.0', 'b': '6356752.3', 'h': '35785863.0', 'lon_0': '140.7', 'proj': 'geos', 'units': 'm'}
Number of columns: 11000
Number of rows: 11000
Area extent: (-5500000.035542117, -5500000.035542117, 5500000.035542117, 5500000.035542117), 'resolution': 1000, 'satellite_altitude': 35783594.40642261, 'calibration': 'reflectance', 'platform_name': 'Himawari-8', 'standard_name': 'toa_bidirectional_reflectance', 'modifiers': (), 'end_time': datetime.datetime(2018, 3, 15, 6, 9, 40, 535246), 'reader': 'ahi_hsd', 'units': '%', 'wavelength': (0.45, 0.47, 0.49), 'satellite_latitude': -0.024062543704014554, 'id': DatasetID(name='B01', wavelength=(0.45, 0.47, 0.49), resolution=1000, polarization=None, calibration='reflectance', modifiers=())}

[djhoese]

Congratulations! You found two bugs.

In your copy of the scmi.py change that ds_info['platform'] to ds_info['platform_name'] and hopefully that should make things work a little better.

When I originally made the SCMI writer for satpy I was using my geocat reader. There is a bug in the geocat reader (typo I guess) where I use platform. Because of that the SCMI writer also uses platform instead of platform_name. Let me know how changing that part of the scmi.py code works for you and we'll go from there.

[WxmanJ]

I was actually just commenting the same thing about "platform" and "platform_name." You stole my thunder. :-P

The writer does work correctly after making the change.

[djhoese]

@WxmanJ Awesome! Thanks for your patience. Does the data load in to AWIPS? It may not look the way you want until style rules are set up, but if it at least loads that's a good sign.

I won't be fixing this bug in the 0.8.x version of satpy, but it will be fixed in future 0.9 releases.

[WxmanJ]

The Wari data does load into AWIPS and I am in the process of setting up style rules.

I will be testing MSG data later today.

@davidh-ssec Kudos on the writer! I am actually quite impressed with the methods you setup for tiling and pixels. That must have taken some time to put together, but it makes life so much easier for us "AWIPS" users in the end. Thank you!

[djhoese]

That must have taken some time to put together

I've probably rewritten it 4 times because of bugs and AWIPS being a pain in what it accepts as input.

Glad it is working now.

[WxmanJ]

David,

I just tried the same type of script using Meteosat-11 data. Upon loading the data into AWIPS, I am getting "No Data" for the display. I have created a style rule for it and it does load. Should I post this as a separate issue?

My traceback below.

[root@satpy ~]# python test_msg.py 
[DEBUG: 2018-03-15 14:41:35 : satpy.scene] Setting 'PPP_CONFIG_DIR' to '/usr/lib/python2.7/site-packages/satpy/etc'
[DEBUG: 2018-03-15 14:41:35 : satpy.readers] Reading ['/usr/lib/python2.7/site-packages/satpy/etc/readers/hrit_msg.yaml', '/usr/lib/python2.7/site-packages/satpy/etc/readers/hrit_msg.yaml', '/usr/lib/python2.7/site-packages/satpy/etc/readers/hrit_msg.yaml']
[DEBUG: 2018-03-15 14:41:36 : satpy.readers.yaml_reader] Assigning to hrit_msg: ['/root/msg_data/H-000-MSG4__-MSG4________-_________-PRO______-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-_________-EPI______-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000003___-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000004___-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000002___-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000007___-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000008___-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000005___-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000001___-201803031200-__', '/root/msg_data/H-000-MSG4__-MSG4________-IR_108___-000006___-201803031200-__']
[INFO: 2018-03-15 14:41:36 : hrit_msg] No IMPF configuration field found in prologue.
[DEBUG: 2018-03-15 14:41:36 : satpy.composites] Looking for composites config file seviri.yaml
[DEBUG: 2018-03-15 14:41:36 : satpy.composites] Looking for composites config file visir.yaml
[DEBUG: 2018-03-15 14:41:36 : hrit_base] Reading time 0:00:00.040485
/usr/lib/python2.7/site-packages/satpy/readers/hrit_msg.py:1000: RuntimeWarning: divide by zero encountered in reciprocal
  data.data[:] **= -1
[DEBUG: 2018-03-15 14:41:36 : hrit_msg] Calibration time 0:00:00.029796
[DEBUG: 2018-03-15 14:41:36 : hrit_base] Reading time 0:00:00.043713
[DEBUG: 2018-03-15 14:41:36 : hrit_msg] Calibration time 0:00:00.042935
[DEBUG: 2018-03-15 14:41:36 : hrit_base] Reading time 0:00:00.036479
[DEBUG: 2018-03-15 14:41:37 : hrit_msg] Calibration time 0:00:00.042631
[DEBUG: 2018-03-15 14:41:37 : hrit_base] Reading time 0:00:00.043688
[DEBUG: 2018-03-15 14:41:37 : hrit_msg] Calibration time 0:00:00.106182
[DEBUG: 2018-03-15 14:41:37 : hrit_base] Reading time 0:00:00.045576
[DEBUG: 2018-03-15 14:41:37 : hrit_msg] Calibration time 0:00:00.057280
[DEBUG: 2018-03-15 14:41:37 : hrit_base] Reading time 0:00:00.038614
[DEBUG: 2018-03-15 14:41:37 : hrit_msg] Calibration time 0:00:00.043964
[DEBUG: 2018-03-15 14:41:37 : hrit_base] Reading time 0:00:00.036673
[DEBUG: 2018-03-15 14:41:37 : hrit_msg] Calibration time 0:00:00.042737
[DEBUG: 2018-03-15 14:41:37 : hrit_base] Reading time 0:00:00.039408
[DEBUG: 2018-03-15 14:41:37 : hrit_msg] Calibration time 0:00:00.035685
{'start_time': datetime.datetime(2018, 3, 3, 12, 0), 'sensor': set(['seviri']), 'end_time': datetime.datetime(2018, 3, 3, 12, 15)}
[INFO: 2018-03-15 14:41:38 : satpy.writers.scmi] Writing product IR_108 to AWIPS SCMI NetCDF file
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] Scaling IR_108 data to fit in netcdf file...
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] Using product valid min 184.219863892 and valid max 331.348510742
[WARNING: 2018-03-15 14:41:38 : satpy.writers.scmi] AWIPS file already exists, will overwrite: AWIPS_AII_MET11_SEVIRI_VIS006_GEOS_01.nc
[INFO: 2018-03-15 14:41:38 : satpy.writers.scmi] Writing tile 'T001' to 'AWIPS_AII_MET11_SEVIRI_VIS006_GEOS_01.nc'
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] Creating dimensions...
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] Creating variables...
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] Creating global attributes...
[WARNING: 2018-03-15 14:41:38 : satpy.writers.scmi] environment ORGANIZATION not set for .production_location attribute, using hostname
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] already have a value for satellite_id
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] Creating projection attributes...
[DEBUG: 2018-03-15 14:41:38 : satpy.writers.scmi] Writing image data...
[INFO: 2018-03-15 14:41:38 : satpy.writers.scmi] writing image data
[DEBUG: 2018-03-15 14:41:39 : satpy.writers.scmi] Writing X/Y navigation data...
[DEBUG: 2018-03-15 14:41:39 : satpy.writers.scmi] y variable shape is (3712,)

[djhoese]

If you do global_scene['VIS006'].show() or whatever dataset you are looking at does it show you what you expect? Or even run global_scene.save_datasets() to save to geotiffs, do they look as you would expect? Meaning, do they have actual data in them?

Edit: Please use 3 back ticks (`) to format your output above.

[WxmanJ]

Yes, data does show. I just loaded an IR108 and saved it as a PNG.

[djhoese]

If you inspect the NetCDF file that is generated by the SCMI writer (I believe it defaults to one tile), what does the data look like? Is it all fills? I have no experience putting this data through the SCMI writer, but it should be have very similarly. Another possibility is that it is not being written with the correct units for either the SCMI writer to convert it to what AWIPS expects or it may be units that AWIPS isn't configured to understand.

[WxmanJ]

The data units are in Kelvin.

Here is a short summary for the IR108:

netcdf AWIPS_AII_MET11_SEVIRI_IR108_GEOS_01 {
dimensions:
	y = 3712 ;
	x = 3712 ;
variables:
	short data(y, x) ;
		data:_FillValue = 32767s ;
		data:coordinates = "y x" ;
		data:scale_factor = 0.004490284f ;
		data:add_offset = 184.2199f ;
		data:units = "kelvin" ;
		data:valid_min = 0s ;
		data:valid_max = 32766s ;
		data:standard_name = "brightness_temperature" ;
		data:grid_mapping = "fixedgrid_projection" ;
	short y(y) ;
		y:scale_factor = -3000.40327858102 ;
		y:add_offset = -5568748.48504637 ;
		y:units = "meters" ;
		y:standard_name = "projection_y_coordinate" ;
	short x(x) ;
		x:scale_factor = -3000.40327858102 ;
		x:add_offset = -5568748.48504637 ;
		x:units = "meters" ;
		x:standard_name = "projection_x_coordinate" ;
	int fixedgrid_projection ;
		fixedgrid_projection:short_name = "some_area_name" ;
		fixedgrid_projection:grid_mapping_name = "geostationary" ;
		fixedgrid_projection:sweep_angle_axis = "y" ;
		fixedgrid_projection:perspective_point_height = 35785831. ;
		fixedgrid_projection:latitude_of_projection_origin = 0.f ;
		fixedgrid_projection:longitude_of_projection_origin = 0.f ;
		fixedgrid_projection:semi_major_axis = 6378169. ;
		fixedgrid_projection:semi_minor_axis = 6356583.8 ;
		fixedgrid_projection:false_easting = 0.f ;
		fixedgrid_projection:false_northing = 0.f ;

// global attributes:
		:_NCProperties = "version=1|netcdflibversion=4.4.1.1|hdf5libversion=1.8.18" ;
		:Conventions = "CF-1.7" ;
		:creator = "UW SSEC - CSPP Polar2Grid" ;
		:creation_time = "2018-03-15T19:15:22" ;
		:physical_element = "IR_108" ;
		:satellite_id = "METEOSAT-11-SEVIRI" ;
		:awips_id = "AWIPS_IR_108" ;
		:sector_id = "Europe" ;
		:tile_row_offset = 0L ;
		:tile_column_offset = 0L ;
		:product_tile_height = 3712L ;
		:product_tile_width = 3712L ;
		:number_product_tiles = 1L ;
		:product_rows = 3712L ;
		:product_columns = 3712L ;
		:pixel_x_size = -3.00040327858102 ;
		:pixel_y_size = -3.00040327858102 ;
		:product_name = "IR_108" ;
		:production_location = "satpy" ;
		:start_date_time = "2018-03-03T12:00:10" ;

If you would like to review the file generated, I am attaching it.
AWIPS_AII_MET11_SEVIRI_IR108_GEOS_01.nc.tar.gz

[djhoese]

I don't have time to look at the file. You'll have to look yourself using either ncdump -v data file.nc or in python. By the looks of the scale factor and offset I would assume the data is fine. If it is all fill values then that explains why AWIPS is seeing no data values everywhere. If some of the values in the NetCDF file are valid then it is an AWIPS issue.

The only other thing that stands out to me is the standard_name which for most other BTs I've seen should be "toa_brightness_temperature". I don't think AWIPS cares about standard_name, but it may be something to look at.

[WxmanJ]

There is some missing data judging by the "_," marks, but for the most part, there is def data there.

[WxmanJ]

Just for the sake of documenting a comparison, here is the print detail from both Himawari-8 and Meteosat-11. Not sure if this will add any value, but I wanted to at least present it. I did notice the DatasetID was not present on Meteosat-11.

Himawari:

ahi/B10: 
	area: On-the-fly area
	calibration: brightness_temperature
	end_time: 2018-03-15 06:09:40.042551
	id: DatasetID(name='B10', wavelength=(7.1, 7.3, 7.5), resolution=2000, polarization=None, calibration='brightness_temperature', modifiers=())
	modifiers: ()
	platform_name: Himawari-8
	polarization: None
	reader: ahi_hsd
	resolution: 2000 m
	satellite_altitude: 35783594.4041
	satellite_latitude: -0.0240619194964
	satellite_longitude: 140.64899483
	standard_name: toa_brightness_temperature
	start_time: 2018-03-15 06:00:20.644159
	units: K
	wavelength: (7.1, 7.3, 7.5) μm
	shape: (5500, 5500)
[[-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 ...
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]]
{'polarization': None, 'sensor': 'ahi', 'start_time': datetime.datetime(2018, 3, 15, 6, 0, 20, 644159), 'name': 'B10', 'satellite_longitude': 140.64899483025076, 'area': Area ID: some_area_name
Description: On-the-fly area
Projection ID: geosh8
Projection: {'a': '6378137.0', 'b': '6356752.3', 'h': '35785863.0', 'lon_0': '140.7', 'proj': 'geos', 'units': 'm'}
Number of columns: 5500
Number of rows: 5500
Area extent: (-5499999.901174725, -5499999.901174725, 5499999.901174725, 5499999.901174725), 'resolution': 2000, 'satellite_altitude': 35783594.404132344, 'calibration': 'brightness_temperature', 'platform_name': 'Himawari-8', 'standard_name': 'toa_brightness_temperature', 'modifiers': (), 'end_time': datetime.datetime(2018, 3, 15, 6, 9, 40, 42551), 'reader': 'ahi_hsd', 'units': 'K', 'wavelength': (7.1, 7.3, 7.5), 'satellite_latitude': -0.024061919496408027, 'id': DatasetID(name='B10', wavelength=(7.1, 7.3, 7.5), resolution=2000, polarization=None, calibration='brightness_temperature', modifiers=())}

Meteosat:

seviri/IR_108: 
	area: On-the-fly area
	calibration: brightness_temperature
	end_time: 2018-03-03 12:12:43.418000
	modifiers: ()
	platform_name: Meteosat-11
	polarization: None
	reader: hrit_msg
	resolution: 3000.40316582 m
	satellite_altitude: 35785831.0
	satellite_latitude: 0.0
	satellite_longitude: 0.0
	standard_name: brightness_temperature
	start_time: 2018-03-03 12:00:10.246000
	units: K
	wavelength: (9.8, 10.8, 11.8) μm
	shape: (3712, 3712)
[[-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 ...
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]]
{'polarization': None, 'modifiers': (), 'name': 'IR_108', 'satellite_longitude': 0.0, 'area': Area ID: some_area_name
Description: On-the-fly area
Projection ID: geosmsg
Projection: {'a': '6378169.0', 'b': '6356583.8', 'h': '35785831.0', 'lon_0': '0.0', 'proj': 'geos', 'units': 'm'}
Number of columns: 3712
Number of rows: 3712
Area extent: (5567248.28340708, 5567248.28340708, -5570248.686685662, -5570248.686685662), 'sensor': 'seviri', 'satellite_altitude': 35785831.0, 'calibration': 'brightness_temperature', 'resolution': 3000.403165817, 'platform_name': 'Meteosat-11', 'standard_name': 'brightness_temperature', 'end_time': datetime.datetime(2018, 3, 3, 12, 12, 43, 418000), 'reader': 'hrit_msg', 'units': 'K', 'wavelength': (9.8, 10.8, 11.8), 'satellite_latitude': 0.0, 'start_time': datetime.datetime(2018, 3, 3, 12, 0, 10, 246000)}

[WxmanJ]

I believe I have resolved the issue with Meteosat-11. I had to run a resample to put it in the correct projection. This may have to do with the inverted raw data.

It now displays in AWIPS.