From 337ef87986f3e9d9eb71b5f1df2c78abd1d8a987 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Thu, 3 Dec 2015 14:45:53 +0100 Subject: [PATCH 1/5] Add FY-3B template config file Signed-off-by: Adam.Dybbroe --- etc/FY-3B.cfg.template | 60 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 etc/FY-3B.cfg.template diff --git a/etc/FY-3B.cfg.template b/etc/FY-3B.cfg.template new file mode 100644 index 00000000..6434e083 --- /dev/null +++ b/etc/FY-3B.cfg.template @@ -0,0 +1,60 @@ +[satellite] +instruments = ('virr',) +satname = fy3 + +[virr-level2] +format=fy3_virr +# clear_FY3B_26293_02-DEC-2015_14-02-13_VIRRX_L1B.HDF +dir=/data/to/my/fengyun3/virr/file/ +filename = clear_FY3B_%(orbit)s_%d-%b-%Y_%H-%M-%S_VIRRX_L1B.HDF + +[virr-1] +name = '1' +frequency = (0.58, 0.63, 0.68) +resolution = 1000 + +[virr-2] +name = '2' +frequency = (0.84, 0.865, 0.89) +resolution = 1000 + +[virr-3] +name = '3' +frequency = (3.55, 3.74, 3.93) +resolution = 1000 + +[virr-4] +name = '4' +frequency = (10.3, 10.8, 11.3) +resolution = 1000 + +[virr-5] +name = '5' +frequency = (11.5, 12.0, 12.5) +resolution = 1000 + +[virr-6] +name = '6' +frequency = (1.55, 1.6, 1.64) +resolution = 1000 + +[virr-7] +name = '7' +frequency = (0.43, 0.455, 0.48) +resolution = 1000 + +[virr-8] +name = '8' +frequency = (0.48, 0.505, 0.53) +resolution = 1000 + +[virr-9] +name = '9' +frequency = (0.53, 0.555, 0.58) +resolution = 1000 + +[virr-10] +name = '10' +frequency = (1.325, 1.36, 1.395) +resolution = 1000 + From 93bfab6fb987b29fea714f3f15919e7b1d54e9be Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Thu, 3 Dec 2015 16:39:14 +0100 Subject: [PATCH 2/5] Apply VIS/NIR calibration including sun-zenith correction Signed-off-by: Adam.Dybbroe --- mpop/satin/fy3_virr.py | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/mpop/satin/fy3_virr.py b/mpop/satin/fy3_virr.py index c563d082..66247d2b 100644 --- a/mpop/satin/fy3_virr.py +++ b/mpop/satin/fy3_virr.py @@ -37,7 +37,7 @@ def load(satscene, *args, **kwargs): """Read data from file and load it into *satscene*. - A possible *calibrate* keyword argument is passed to the AAPP reader. + A possible *calibrate* keyword argument is passed to the AAPP reader. Should be 0 for off (counts), 1 for default (brightness temperatures and reflectances), and 2 for radiances only. @@ -94,8 +94,18 @@ def load_virr(satscene, options): calibrate = options['calibrate'] - # Get the calibration information h5f = h5py.File(filename, 'r') + + # Get geolocation information + lons = h5f['Longitude'][:] + lats = h5f['Latitude'][:] + sunz = h5f['SolarZenith'][:] + slope = h5f['SolarZenith'].attrs['Slope'][0] + intercept = h5f['SolarZenith'].attrs['Intercept'][0] + sunz = sunz * slope + intercept + sunz = np.where(np.greater(sunz, 85.0), 85.0, sunz) + + # Get the calibration information # Emissive radiance coefficients: emis_offs = h5f['Emissive_Radiance_Offsets'][:] emis_scales = h5f['Emissive_Radiance_Scales'][:] @@ -103,13 +113,16 @@ def load_virr(satscene, options): # VIS/NIR calibration stuff: refsb_cal_coeff = h5f.attrs['RefSB_Cal_Coefficients'] + visnir_scales = refsb_cal_coeff[0::2] + visnir_offs = refsb_cal_coeff[1::2] + refsb_effective_wl = h5f.attrs['RefSB_Effective_Wavelength'] + # Read the band data: for dset in datasets: band_data = h5f[dset] valid_range = band_data.attrs['valid_range'] LOGGER.debug("valid-range = " + str(valid_range)) - # FIXME! There seem to be useful data outside the valid range! fillvalue = band_data.attrs['_FillValue'] band_names = band_data.attrs['band_name'].split(',') slope = band_data.attrs['Slope'] @@ -129,18 +142,18 @@ def load_virr(satscene, options): bandmask = np.logical_or(np.less(data, valid_range[0]), np.greater(data, valid_range[1])) - # if calibrate: - # data = slopes[mersi_band_index] * ( - # data - np.array([sv_dn_average[mersi_band_index]]).transpose()) + if calibrate: + if dset in ['EV_Emissive']: + data = (np.array([emis_offs[:, i]]).transpose() + + data * np.array([emis_scales[:, i]]).transpose()) + if dset in ['EV_RefSB']: + data = (visnir_offs[i] + + data * visnir_scales[i]) / np.cos(np.deg2rad(sunz)) satscene[band] = np.ma.masked_array(data, mask=bandmask, copy=False) - # Get geolocation information - lons = h5f['Longitude'][:] - lats = h5f['Latitude'][:] - pdb.set_trace() from pyresample import geometry satscene.area = geometry.SwathDefinition(lons=lons, lats=lats) From 4568f894c517b7299da6e095db263939cd5bec9b Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Thu, 3 Dec 2015 16:40:26 +0100 Subject: [PATCH 3/5] Remove old FY3 mersi reader Signed-off-by: Adam.Dybbroe --- mpop/satin/fy3_mersi_aggr1km.py | 162 -------------------------------- 1 file changed, 162 deletions(-) delete mode 100644 mpop/satin/fy3_mersi_aggr1km.py diff --git a/mpop/satin/fy3_mersi_aggr1km.py b/mpop/satin/fy3_mersi_aggr1km.py deleted file mode 100644 index cadb222c..00000000 --- a/mpop/satin/fy3_mersi_aggr1km.py +++ /dev/null @@ -1,162 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# Copyright (c) 2010. - -# SMHI, -# Folkborgsvägen 1, -# Norrköping, -# Sweden - -# Author(s): - -# Adam Dybbroe - -# This file is part of mpop. - -# mpop 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. - -# mpop 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 -# mpop. If not, see . - -"""Interface to 1km agregated MERSI hdf5 level 1 format. -""" -import os.path -from ConfigParser import ConfigParser - -import _pyhl -import numpy as np - -from mpop import BASE_PATH - - -def load(satscene): - """Read data from file and load it into *satscene*. - """ - conf = ConfigParser() - conf.read(os.path.join(BASE_PATH, "etc", satscene.fullname + ".cfg")) - options = {} - for option, value in conf.items(satscene.instrument_name+'-level2', raw = True): - options[option] = value - CASES[satscene.instrument_name](satscene, options) - - -def load_1km_aggregated_mersi(satscene, options): - """Read 1km agregated mersi data from file and load it into *satscene*. - """ - # Example: FY3A_MERSI_GBAL_L1_20100308_0915_1000M_MS.HDF - filename = satscene.time_slot.strftime("FY3A_MERSI_GBAL_L1_%Y%m%d_%H%M_1000M_MS.HDF") - filename = os.path.join(options["dir"], filename) - - a = _pyhl.read_nodelist(filename) - b = a.getNodeNames() - # Should only select/fetch the datasets needed. FIXME! - a.selectAll() - a.fetch() - - # MERSI Channel 1-4: EV_250_Aggr.1KM_RefSB - # MERSI Channel 5: EV_250_Aggr.1KM_Emissive - # MERSI Channel 6-20: EV_1KM_RefSB - - datasets = ['/EV_250_Aggr.1KM_RefSB', - '/EV_250_Aggr.1KM_Emissive', - '/EV_1KM_RefSB'] - - for nodename in datasets: - band_data = a.getNode(nodename).data() - valid_range = a.getNode('%s/valid_range' % (nodename)).data() - band_names = a.getNode('%s/band_name' % (nodename)).data().split(",") - # Special treatment of the situation where the bandnames are stored - # as '6~20' (quite inconvenient): - if '6~20' in band_names: - band_names = ['6','7','8','9','10','11','12','13', - '14','15','16','17','18','19','20'] - - for (i, band) in enumerate(band_names): - if band not in satscene.channels_to_load: - continue - - satscene[band] = np.ma.masked_outside(band_data[i], - valid_range[0], - valid_range[1], - copy = False) - satscene[band].info = { - 'var_name' : 'ch'+str(band), - 'var_data' : satscene[band].data, - 'var_dim_names': ('x','y'), - '_FillValue' : 32767, - 'standard_name' : '', - 'short_name' : band, - 'scale_factor' : 1.0, - 'add_offset' : 0.0, - } - - - satscene.info = { - 'var_children' : [ #{'var_name' : 'lat', 'var_callback': Functor(satscene.get_lat, low_res), 'var_dim_names': ('x','y') }, - #{'var_name' : 'lon', 'var_callback' : Functor(satscene.get_lon, low_res) , 'var_dim_names': ('x','y') }, - ## {'var_name' : 'lat_hrvis', 'var_data' : satscene.lat[high_res]}, - ## {'var_name' : 'lon_hrvis', 'var_data' : satscene.lon[high_res]}, - ], - 'Satellite' : satscene.fullname, - 'Antenna' : 'None', - 'Receiver' : 'SMHI' , - 'Time' : satscene.time_slot.strftime("%Y-%m-%d %H:%M:%S UTC"), - #'Area_Name' : satscene.area_id, - 'Area_Name' : "swath", - 'Projection' : 'satproj', - 'Platform' : 'fengyun', - 'Number' : '3a', - 'Service' : '', - #'Columns' : satscene.channels[0].shape[1], - #'Lines' : satscene.channels[0].shape[0], - 'SampleX' : 1.0, - 'SampleY' : 1.0, - 'title' : 'MERSI Level 1', - } - - -def get_lat_lon(satscene, resolution): - """Read lat and lon. - """ - del resolution - - conf = ConfigParser() - conf.read(os.path.join(BASE_PATH, "etc", satscene.fullname + ".cfg")) - options = {} - for option, value in conf.items(satscene.instrument_name+'-level2', raw = True): - options[option] = value - - return LAT_LON_CASES[satscene.instrument_name](satscene, options) - -def get_lat_lon_1km_aggregated_mersi(satscene, options): - """Read latitude and longitude for each (aggregated) pixel. - """ - # Example: FY3A_MERSI_GBAL_L1_20100308_0915_1000M_MS.HDF - filename = satscene.time_slot.strftime("FY3A_MERSI_GBAL_L1_%Y%M%D_%H%M_1000M_MS.HDF") - filename = os.path.join(options["dir"], filename) - - a = _pyhl.read_nodelist(filename) - b = a.getNodeNames() - # Should only select/fetch the datasets needed. FIXME! - a.selectAll() - a.fetch() - - lat = a.getNode("/Latitude").data() - lon = a.getNode("/Longitude").data() - - return lat, lon - -CASES = { - "mersi": load_1km_aggregated_mersi - } - -LAT_LON_CASES = { - "mersi": get_lat_lon_1km_aggregated_mersi - } From 3fc9ab87d249440ae39a1e225ce663eb1eedec32 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Fri, 4 Dec 2015 10:47:08 +0100 Subject: [PATCH 4/5] Update EARS config files for new (2014) PPS product format Signed-off-by: Adam.Dybbroe --- etc/EARSMetop-A.cfg.template | 6 ++-- etc/EARSMetop-B.cfg.template | 62 ++++++++++++++++++++++++++++++++++++ etc/EARSNOAA-19.cfg.template | 8 ++--- 3 files changed, 69 insertions(+), 7 deletions(-) create mode 100644 etc/EARSMetop-B.cfg.template diff --git a/etc/EARSMetop-A.cfg.template b/etc/EARSMetop-A.cfg.template index 25e10dcd..089fd8b8 100644 --- a/etc/EARSMetop-A.cfg.template +++ b/etc/EARSMetop-A.cfg.template @@ -3,9 +3,9 @@ variant = EARS instruments = ('avhrr/3',) [avhrr/3-level3] -filename = %(product)s_%Y%m%d_%H%M00_%(satellite)s.h5* -dir = /path/to/my/ears/nwc -format = nwcsaf_pps +format = nc_pps_l2.PPSReader +cloud_product_filename = W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,METOPA+%(product)s_C_EUMS_%Y%m%d%H%M00_%(orbit)s.nc* +cloud_product_geofilename = W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,METOPA+CMA_C_EUMS_%Y%m%d%H%M00_%(orbit)s.nc* [avhrr/3-level1] format=eps1a diff --git a/etc/EARSMetop-B.cfg.template b/etc/EARSMetop-B.cfg.template new file mode 100644 index 00000000..127b2993 --- /dev/null +++ b/etc/EARSMetop-B.cfg.template @@ -0,0 +1,62 @@ +[satellite] +satname = Metop-B +number = +variant = EARS +instruments = ('avhrr/3',) + +[avhrr/3-level3] +format = nc_pps_l2.PPSReader +# W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,NOAA19+CTTH_C_EUMS_20150820141700_33658.nc.bz2 +cloud_product_filename = W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,METOPB+%(product)s_C_EUMS_%Y%m%d%H%M00_%(orbit)s.nc* +cloud_product_geofilename = W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,METOPB+CMA_C_EUMS_%Y%m%d%H%M00_%(orbit)s.nc* + +[avhrr/3-level1] +format=eps1a +shortname=M01 +dir=/path/to/my/ears/avhrr +filename=AVHR_HRP_00_M01_%Y%m%d%H%M* + +[avhrr/3-granules] +type=bzipped eps1a +granularity=60 +full_scan_period=0.1667 +scan_width=2048 +dir=/path/to/my/ears/avhrr +filename=AVHR_HRP_00_M02_%Y%m%d%H%M* + + +[avhrr/3-1] +name = '1' +frequency = (0.58, 0.63, 0.68) +resolution = 1090 +size = (2048,) + +[avhrr/3-2] +name = '2' +frequency = (0.725, 0.8625, 1.0) +resolution = 1090 +size = (2048,) + +[avhrr/3-3] +name = '3A' +frequency = (1.58, 1.61, 1.64) +resolution = 1090 +size = (2048,) + +[avhrr/3-4] +name = '3B' +frequency = (3.55, 3.74, 3.93) +resolution = 1090 +size = (2048,) + +[avhrr/3-5] +name = '4' +frequency = (10.3, 10.8, 11.3) +resolution = 1090 +size = (2500, 2500) + +[avhrr/3-6] +name = '5' +frequency = (11.5, 12.0, 12.5) +resolution = 1090 +size = (2500, 2500) diff --git a/etc/EARSNOAA-19.cfg.template b/etc/EARSNOAA-19.cfg.template index 8e566b4d..e0dcc5b5 100644 --- a/etc/EARSNOAA-19.cfg.template +++ b/etc/EARSNOAA-19.cfg.template @@ -3,9 +3,10 @@ variant = EARS instruments = ('avhrr/3',) [avhrr/3-level3] -filename = %(product)s_%Y%m%d_%H%M00_%(satellite)s.h5* -dir = /path/to/my/ears/nwc -format = nwcsaf_pps +format = nc_pps_l2.PPSReader +# W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,NOAA19+CTTH_C_EUMS_20150820141700_33658.nc.bz2 +cloud_product_filename = W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,NOAA19+%(product)s_C_EUMS_%Y%m%d%H%M00_%(orbit)s.nc* +cloud_product_geofilename = W_XX-EUMETSAT-Darmstadt,SING+LEV+SAT,NOAA19+CMA_C_EUMS_%Y%m%d%H%M00_%(orbit)s.nc* [avhrr/3-level1] format=hrpt @@ -21,7 +22,6 @@ scan_width=2048 dir=/path/to/my/ears/avhrr filename=avhrr_%Y%m%d_%H%M%S_noaa19.hrp.bz2 - [avhrr/3-1] name = '1' frequency = (0.58, 0.63, 0.68) From c4168f1048e1d1dbeebb14a75ca0fae34906907c Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Fri, 4 Dec 2015 15:18:38 +0100 Subject: [PATCH 5/5] Add brightness temperature calibration to the IR bands Signed-off-by: Adam.Dybbroe --- mpop/satin/fy3_virr.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/mpop/satin/fy3_virr.py b/mpop/satin/fy3_virr.py index 66247d2b..c673e6fb 100644 --- a/mpop/satin/fy3_virr.py +++ b/mpop/satin/fy3_virr.py @@ -30,7 +30,7 @@ from ConfigParser import ConfigParser from mpop import CONFIG_PATH import h5py -import pdb +from pyspectral.blackbody import blackbody_wn_rad2temp as rad2temp LOGGER = logging.getLogger('virr') @@ -109,6 +109,10 @@ def load_virr(satscene, options): # Emissive radiance coefficients: emis_offs = h5f['Emissive_Radiance_Offsets'][:] emis_scales = h5f['Emissive_Radiance_Scales'][:] + + # Central wave number (unit = cm-1) for the three IR bands + # It is ordered according to decreasing wave number (increasing wavelength): + # 3.7 micron, 10.8 micron, 12 micron emiss_centroid_wn = h5f.attrs['Emmisive_Centroid_Wave_Number'] # VIS/NIR calibration stuff: @@ -146,6 +150,12 @@ def load_virr(satscene, options): if dset in ['EV_Emissive']: data = (np.array([emis_offs[:, i]]).transpose() + data * np.array([emis_scales[:, i]]).transpose()) + # Radiance to Tb conversion. + # Pyspectral wants SI units, + # but radiance data are in mW/m^2/str/cm^-1 and wavenumbers are in cm^-1 + # Therefore multply wavenumber by 100 and radiances by + # 10^-5 + data = rad2temp(emiss_centroid_wn[i] * 100., data * 1e-5) if dset in ['EV_RefSB']: data = (visnir_offs[i] + data * visnir_scales[i]) / np.cos(np.deg2rad(sunz))