diff --git a/etc/hsaf10.cfg.template b/etc/hsaf10.cfg.template new file mode 100755 index 00000000..191bdd08 --- /dev/null +++ b/etc/hsaf10.cfg.template @@ -0,0 +1,26 @@ +[satellite] +satname = 'hsaf' +projection = 'geos(0.0)' +number = '10' +instruments = ('seviri',) +proj4_params = 'proj=geos lon_0=0.0 lat_0=0.00 a=6378169.00 b=6356583.80 h=35785831.00' + +#[seviri-level1] +##filename = h03_%Y%m%d_%H%M_rom.grb.nc4 +#filename = h03_%Y%m%d_%H%M_rom.* +#format = 'read_h03' +##dir = /data/COALITION2/database/meteosat/HSAF/%Y/%m/%d +#dir = /data/cinesat/in/hsaf + +[seviri-level2] +#filename = h03_%Y%m%d_%H%M_rom.grb.nc4 +filename = h03_%Y%m%d_%H%M_rom.* +format = 'hsaf_h03' +#dir = /data/COALITION2/database/meteosat/HSAF/%Y/%m/%d +dir = /data/cinesat/in/hsaf + +[seviri-1] +frequency = (0.00, 0.00, 0.00) +resolution = 3000.403165817 +name = 'h03' +size = (1900, 900) \ No newline at end of file diff --git a/etc/odyssey.cfg.template b/etc/odyssey.cfg.template new file mode 100755 index 00000000..cf7c9653 --- /dev/null +++ b/etc/odyssey.cfg.template @@ -0,0 +1,25 @@ +[satellite] +satname = 'odyssey' +variant = +number = +instruments = ('radar',) + +[radar-level2] +dir = /data/cinesat/in/radar +format = odyssey_radar +projection = odyssey + +[RATE] +#filename = T_PAAH21_C_EUOC_%Y%m%d%H%M??.hdf +filename = meteoswiss.radar.euoc_rain_rate.%Y%m%d%H%M.hdf5 +name = 'RATE' + +[DBZH] +#filename = T_PABH21_C_EUOC_%Y%m%d%H%M??.hdf +filename = meteoswiss.radar.euoc_maximum_reflectivit.%Y%m%d%H%M.hdf5 +name = 'DBZH' + +[ACRR] +#filename = T_PASH21_C_EUOC_%Y%m%d%H%M??.hdf +filename = meteoswiss.radar.euoc_1h_accumulation.%Y%m%d%H%M.hdf5 +name = 'ACRR' diff --git a/mpop/satin/hsaf_h03.py b/mpop/satin/hsaf_h03.py new file mode 100644 index 00000000..7243a5bc --- /dev/null +++ b/mpop/satin/hsaf_h03.py @@ -0,0 +1,226 @@ +""" +Reader for EUMETSATs Hydrology SAF (HSAF) h03 product +HSAF website http://hsaf.meteoam.it +h03 product is precipitation rate at the ground +by GEO(MSG)/Infrared supported by LEO/Microwave +http://hsaf.meteoam.it/precipitation.php?tab=3 + +After registration the data is available from +ftp://ftphsaf.meteoam.it/h03 + +possible accepted formats for this reader are: +* grib as provided by HSAF +* netCDF (grib file converted with cdo) + +- Initial version: + 2015-07-23 Ulrich Hamann (MeteoSwiss) +""" + +from ConfigParser import ConfigParser +from mpop import CONFIG_PATH +import os +import numpy.ma as ma +from glob import glob +import datetime + +def load(satscene, **kargs): + """Reader for EUMETSATs Hydrology SAF (HSAF) h03 product + h03 product is precipitation rate at the ground + by GEO(MSG)/Infrared supported by LEO/Microwave + http://hsaf.meteoam.it/precipitation.php?tab=3 + """ + + # Read config file content + conf = ConfigParser() + conf.read(os.path.join(CONFIG_PATH, satscene.fullname + ".cfg")) + values = {"orbit": satscene.orbit, + "satname": satscene.satname, + "number": satscene.number, + "instrument": satscene.instrument_name, + "satellite": satscene.fullname + } + + # end of scan time 12min after start + end_time = satscene.time_slot + datetime.timedelta(minutes=12) + + filepath = end_time.strftime(conf.get("seviri-level2", "dir",raw=True)) + filepattern = end_time.strftime(conf.get("seviri-level2", "filename",raw=True)) % values + filename = os.path.join( filepath, filepattern) + + print "... search for file: ", filename + filenames=glob(str(filename)) + if len(filenames) == 0: + print "*** Error, no file found" + quit() + elif len(filenames) > 1: + print "*** Warning, more than 1 datafile found: " + for filename in filenames: + print " ", filename + + # possible formats: h03_20150513_1557_rom.grb.gz, h03_20150513_1612_rom.grb, h03_20150513_1612_rom.nc + fileformats = [filename.split(".")[-1] for filename in filenames] + # try to find grb file + + if 'grb' in fileformats: + # read grib + data, fill_value, units, long_name = read_h03_grib(filenames[fileformats.index('grb')]) + elif 'nc' in fileformats: + # read netCDF + data, fill_value, units, long_name = read_h03_netCDF(filenames[fileformats.index('nc')]) + elif 'gz' in fileformats: + # unzip + from subprocess import call + infile = filenames[fileformats.index('gz')] + outfile = infile[:-3] + print " unizp ", infile + # gunzip -c h03_20150513_1557_rom.grb.gz > h03_20150513_1557_rom.grb + # call("/bin/gunzip "+ infile +" 2>&1", shell=True) # dont keep gz file + call("/bin/gunzip -c "+ infile+" > "+ outfile +" 2>&1", shell=True) # keep gz file + # check format of gunziped file + if outfile.split(".")[-1] == 'grb': + data, fill_value, units, long_name = read_h03_grib(outfile) + elif outfile.split(".")[-1] == 'nc': + data, fill_value, units, long_name = read_h03_netCDF(outfile) + + if units == "kg m**-2 s**-1" or units == "kg m-2s-1": + data *= 3600 + units = "kg m-2 h-1" + + satscene['h03'] = data + satscene['h03'].fill_value = fill_value + satscene['h03'].units = units + satscene['h03'].long_name = long_name + satscene['h03'].product_name = 'h03' + + # personal communication with help desk + # Each H03 grib file contains precipitation data of a 900x1900 pixel sub-area of the SEVIRI full disk area (3712x3712 pixels). + # The first pixel of H03 (pixel (1,1)) grib file corresponds to Seviri pixel (1095,85) if the Seviri pixel (1,1) is in the Nort-East. + # I can confirm that only the prime satellite is used (position subsatellite longitude 0 degree East). + # For the future we are thinking to disseminate the h03 outputs already corrected in parallax. + + # conversion of above information to correct AreaDefinition + # full_disk = get_area_def("SeviriDiskFull") + # from mpop.projector import get_area_def + # import numpy as np + # np.array(area_def.get_proj_coords(data_slice=(85+900,1095 ))) - 3000.40316582 / 2. + # array([-2284807.01076965, 2611850.9558437 ]) + # np.array(area_def.get_proj_coords(data_slice=(85 ,1095+1900))) + 3000.40316582 / 2. + # array([ 3418959.40744847, 5315214.20824482]) + # or + # aex = full_disk.get_area_extent_for_subsets(985,1095,85,2995) + + proj = {'proj': 'geos', 'a': '6378169.0', 'b': '6356583.8', 'h': '35785831.0', 'lon_0': '0.0'} + aex = (-2284807.01076965, 2611850.9558437, 3418959.40744847, 5315214.20824482) + + from pyresample.geometry import AreaDefinition + satscene.area = AreaDefinition("hsaf", + "hsaf", + "geos0", + proj, + 1900, + 900, + aex) + +def read_h03_grib(filename): + + try: + import pygrib + except ImportError: + print "... module pygrib needs to be installed" + quit() + + # see http://pygrib.googlecode.com/svn/trunk/docs/pygrib-module.html + print("... read data from %s" % str(filename)) + + grbs = pygrib.open(filename) + + #print(grbs) + + #print 'inventory' + #for grb in grbs: + # print(grb) + #print 'end inventory' + + long_name = 'Instantaneous rain rate' + units = 'kg m**-2 s**-1' + _FillValue = 0.0 + grb = grbs.select(name=long_name)[0] + # print(grb) + data = ma.asarray(grb.values) + data.mask = (data == 0.0) + + print ' fill_value: ', 0 + print ' units: ', units + print ' long_name: ', long_name + print ' datatype: ', type(data) + print ' shape: ', data.shape + print ' min/max: ', data.min(), data.max() + + return data, _FillValue, units, long_name + + +def read_h03_netCDF(filename): + + try: + from netCDF4 import Dataset + except ImportError: + print "... module netCDF4 needs to be installed" + quit() + + print("... read data from %s" % str(filename)) + + # Load data from netCDF file + # see also http://netcdf4-python.googlecode.com/svn/trunk/docs/netCDF4-module.html + ds = Dataset(filename, 'r') + + if 'irrate' in ds.variables: + print ' found variable irrate' + var_name='irrate' # converted with: cdo -f nc4c -z zip copy infile outfile + elif 'IRRATE_P30_GSV0' in ds.variables: + print ' found variable IRRATE_P30_GSV0' + var_name='IRRATE_P30_GSV0' # converted with: ncl_convert2nc h03_20150529_0827_rom.grb -e grb -nc4c -cl 9 + #variables: + # float IRRATE_P30_GSV0(ygrid_0, xgrid_0) ; + # IRRATE_P30_GSV0:initial_time = "05/29/2015 (08:27)" ; + # IRRATE_P30_GSV0:parameter_template_discipline_category_number = 30, 3, 1, 1 ; + # IRRATE_P30_GSV0:parameter_discipline_and_category = "Space products, Quantitative products" ; + # IRRATE_P30_GSV0:grid_type = "Space view perspective or orthographic" ; + # IRRATE_P30_GSV0:_FillValue = 1.e+20f ; + # IRRATE_P30_GSV0:units = "kg m-2s-1" ; + # IRRATE_P30_GSV0:long_name = "Instantaneous rain rate" ; + # IRRATE_P30_GSV0:production_status = "Operational test products" ; + # IRRATE_P30_GSV0:center = "Rome (RSMC)" ; + print '*** Error, does not work for unknown reason' + print ' data.mask = (data == _FillValue) | (data == 0.0) produce error' + quit() + + #print type(ds.variables[var_name]) + #print dir(ds.variables[var_name]) + + _FillValue = ds.variables[var_name]._FillValue + # or fill_value = ds.variables[var_name].getncattr('_FillValue') + units = ds.variables[var_name].units + long_name = ds.variables[var_name].long_name + + # Read variable corresponding to channel name + data = ma.asarray(ds.variables[var_name]) + + print ' fill_value: ', ds.variables[var_name]._FillValue + print ' units: ', ds.variables[var_name].units + print ' long_name: ', ds.variables[var_name].long_name + print ' datatype: ', ds.variables[var_name].datatype + print ' shape: ', data.shape + print ' min/max: ', data.min(), data.max() + + if len(data.shape) == 3: + if data.shape[0] == 1: + print " reduce to 2 dimensions (skip time dimension)" + data = ma.asarray(ds.variables[var_name][0,:,:]) + else: + print "*** Error, unknown netCDF file format in h03_nc.py" + print " probably more time steps in one file (not implemented yet)" + quit() + + data.mask = (data == _FillValue) | (data == 0.0) + + return data, _FillValue, units, long_name diff --git a/mpop/satin/msg_seviri_hdf.py b/mpop/satin/msg_seviri_hdf.py index 8b99cfd0..60e18ea7 100755 --- a/mpop/satin/msg_seviri_hdf.py +++ b/mpop/satin/msg_seviri_hdf.py @@ -79,10 +79,10 @@ def load(satscene, calibrate=True, area_extent=None, **kwargs): dt_end = 12 else: from my_msg_module import check_RSS - RSS = check_RSS(satscene.number, satscene.time_slot) + RSS = check_RSS(satscene.sat_nr(), satscene.time_slot) if RSS == None: print "*** Error in mpop/satin/msg_seviri_hdf.py" - print " satellite MSG", satscene.number ," is not active yet" + print " satellite MSG", satscene.sat_nr() ," is not active yet" quit() else: if RSS: @@ -206,7 +206,7 @@ def load(satscene, calibrate=True, area_extent=None, **kwargs): # satellite ID number hdr["SatelliteDefinition"] = dict() - hdr["SatelliteDefinition"]["SatelliteId"] = SatelliteIds[satscene.number] + hdr["SatelliteDefinition"]["SatelliteId"] = SatelliteIds[str(satscene.sat_nr())] # processing hdr["Level 1_5 ImageProduction"] = dict() diff --git a/mpop/satin/odyssey_radar.py b/mpop/satin/odyssey_radar.py new file mode 100644 index 00000000..d71775bf --- /dev/null +++ b/mpop/satin/odyssey_radar.py @@ -0,0 +1,222 @@ +import Image +import glob +import os +from ConfigParser import ConfigParser +import numpy as np +import numpy.ma as ma +from mpop import CONFIG_PATH + +import pyresample +import logging + +import h5py + +LOG = logging.getLogger(__name__) + +ODIM_H5_FIELD_NAMES = { + 'TH': 'total_power', # uncorrected reflectivity, horizontal + 'TV': 'total_power', # uncorrected reflectivity, vertical + 'DBZH': 'reflectivity', # corrected reflectivity, horizontal + 'DBZV': 'reflectivity', # corrected reflectivity, vertical + 'ZDR': 'differential_reflectivity', # differential reflectivity + 'RHOHV': 'cross_correlation_ratio', + 'LDR': 'linear_polarization_ratio', + 'PHIDP': 'differential_phase', + 'KDP': 'specific_differential_phase', + 'SQI': 'normalized_coherent_power', + 'SNR': 'signal_to_noise_ratio', + 'VRAD': 'velocity', + 'WRAD': 'spectrum_width', + 'QIND': 'quality_index', + 'RATE': 'precip', # precip + 'ACRR': 'accu_precip', # 1 hour ACCU +} + + +def load(satscene, *args, **kwargs): + """Loads the *channels* into the satellite *scene*. + """ + # + # Dataset information + # + # Read config file content + conf = ConfigParser() + conf.read(os.path.join(CONFIG_PATH, satscene.fullname + ".cfg")) + + values = {"orbit": satscene.orbit, + "satname": satscene.satname, + "number": satscene.number, + "instrument": satscene.instrument_name, + "satellite": satscene.fullname + } + + # projection info + projectionName = conf.get("radar-level2", "projection") + projection = pyresample.utils.load_area(os.path.join(CONFIG_PATH, "areas.def"), projectionName) + satscene.area = projection + + for chn_name in satscene.channels_to_load: + filename = os.path.join( + satscene.time_slot.strftime(conf.get("radar-level2", "dir", raw=True)) % values, + satscene.time_slot.strftime(conf.get(chn_name, "filename", raw=True)) % values ) + + # Load data from the h5 file + LOG.debug("filename: "+filename) + filenames=glob.glob(str(filename)) + + if len(filenames) == 0: + LOG.debug("no input file found: "+filename) + print "no input file found:"+filename + quit() + else: + filename = glob.glob(str(filename))[0] + + # open the file + hfile = h5py.File(filename, 'r') + odim_object = hfile['what'].attrs['object'] + if odim_object != 'COMP': + raise NotImplementedError('object: %s not implemented.' % (odim_object)) + else: + # File structure + + #>>> hfile.keys() + #[u'dataset1', u'dataset2', u'how', u'what', u'where'] + + + #>>> for f in hfile['what'].attrs.keys(): + #... print "hfile['what'].attrs['",f,"']=", hfile['what'].attrs[f] + # + #hfile['what'].attrs[' object ']= COMP + #hfile['what'].attrs[' version ']= H5rad 2.0 + #hfile['what'].attrs[' date ']= 20151201 + #hfile['what'].attrs[' time ']= 060000 + #hfile['what'].attrs[' source ']= ORG:247 + + #>>> for f in hfile['where'].attrs.keys(): + #... print "hfile['where'].attrs['",f,"']=", hfile['where'].attrs[f] + # + #hfile['where'].attrs[' projdef ']= +proj=laea +lat_0=55.0 +lon_0=10.0 +x_0=1950000.0 +y_0=-2100000.0 +units=m +ellps=WGS84 + #hfile['where'].attrs[' xsize ']= 1900 + #hfile['where'].attrs[' ysize ']= 2200 + #hfile['where'].attrs[' xscale ']= 2000.0 + #hfile['where'].attrs[' yscale ']= 2000.0 + #hfile['where'].attrs[' LL_lon ']= -10.4345768386 + #hfile['where'].attrs[' LL_lat ']= 31.7462153193 + #hfile['where'].attrs[' UL_lon ']= -39.5357864125 + #hfile['where'].attrs[' UL_lat ']= 67.0228327583 + #hfile['where'].attrs[' UR_lon ']= 57.8119647501 + #hfile['where'].attrs[' UR_lat ']= 67.6210371028 + #hfile['where'].attrs[' LR_lon ']= 29.4210386356 + #hfile['where'].attrs[' LR_lat ']= 31.9876502779 + + # hfile['how'].attrs['nodes'] + # list of radar in composite + + #>>> for f in hfile['dataset1']['what'].attrs.keys(): + #... print "hfile['dataset1'][what].attrs['",f,"']=", hfile['dataset1']['what'].attrs[f] + # + #hfile['dataset1'][what].attrs[' product ']= COMP + #hfile['dataset1'][what].attrs[' startdate ']= 20151201 + #hfile['dataset1'][what].attrs[' starttime ']= 055000 + #hfile['dataset1'][what].attrs[' enddate ']= 20151201 + #hfile['dataset1'][what].attrs[' endtime ']= 060500 + #hfile['dataset1'][what].attrs[' quantity ']= RATE + #hfile['dataset1'][what].attrs[' gain ']= 1.0 + #hfile['dataset1'][what].attrs[' offset ']= 0.0 + #hfile['dataset1'][what].attrs[' nodata ']= -9999000.0 + #hfile['dataset1'][what].attrs[' undetect ']= -8888000.0 + #>>> for f in hfile['dataset2']['what'].attrs.keys(): + #... print "hfile['dataset2'][what].attrs['",f,"']=", hfile['dataset2']['what'].attrs[f] + # + #hfile['dataset2'][what].attrs[' product ']= COMP + #hfile['dataset2'][what].attrs[' startdate ']= 20151201 + #hfile['dataset2'][what].attrs[' starttime ']= 055000 + #hfile['dataset2'][what].attrs[' enddate ']= 20151201 + #hfile['dataset2'][what].attrs[' endtime ']= 060500 + #hfile['dataset2'][what].attrs[' quantity ']= QIND + #hfile['dataset2'][what].attrs[' gain ']= 1.0 + #hfile['dataset2'][what].attrs[' offset ']= 0.0 + #hfile['dataset2'][what].attrs[' nodata ']= -9999000.0 + #hfile['dataset2'][what].attrs[' undetect ']= -8888000.0 + + _xsize = hfile['where'].attrs['xsize'] + _ysize = hfile['where'].attrs['ysize'] + #nbins= _xsize * _ysize + + #projection = hfile['where'].attrs['projdef'] + + datasets = [k for k in hfile if k.startswith('dataset')] + datasets.sort() + nsweeps = len(datasets) + + try: + ds1_what = hfile[datasets[0]]['what'].attrs + except KeyError: + # if no how group exists mock it with an empty dictionary + ds1_what = {} + + _type = '' + if 'product' in ds1_what: + LOG.debug("product: "+ds1_what['product']) + if ds1_what['product'] == 'COMP': + if 'quantity' in ds1_what: + _type = ds1_what['quantity'] + LOG.debug("product_type: "+_type) + + #for chn_name in satscene.channels_to_load: + # if chn_name == _type: + + raw_data = hfile[datasets[0]]['data1']['data'][:] + raw_data = raw_data.reshape(_ysize,_xsize) + + # flag no data + if 'nodata' in ds1_what: + nodata = ds1_what['nodata'] + data = np.ma.masked_equal(raw_data, nodata) + else: + data = np.ma.masked_array(raw_data) + + mask = np.ma.masked_array( raw_data == nodata ) + mask = np.ma.masked_equal( mask, False) + + # flag undetect data + if 'undetect' in ds1_what: + undetect = ds1_what['undetect'] + data[data == undetect] = np.ma.masked + + #from trollimage.image import Image as trollimage + #img = trollimage(mask, mode="L", fill_value=[1,1,1]) # [0,0,0] [1,1,1] + #from trollimage.colormap import rainbow + #img.colorize(rainbow) + #img.show() + #quit() + + # gain/offset adjustment + if 'offset' in ds1_what: + offset = ds1_what['offset'] + else: + offset = 0.0 + + if 'gain' in ds1_what: + gain = ds1_what['gain'] + else: + gain = 1.0 + + data *= gain + offset + + satscene[chn_name] = data + satscene[chn_name+'-MASK'] = mask + + LOG.debug(" *** channel:"+chn_name) + + if _type == 'DBZH': + units = 'dBZ' + + if _type == 'RATE': + units = 'mm/h' + + if _type == 'ACRR': + units = 'mm' + + satscene[chn_name].info["units"] = units + LOG.debug("channel:"+chn_name+" units:"+units)