In [1]:
#PROBLEM: reflectance estimated with demo_refl039 & get_reflectance is a bit different
#sza (from pytroll) & sunz (from knmi) is quite a bit different
#so far missing: cldmask, cph, qa, structural params...
In [1]:
import datetime
import logging
from mpop.satellites import GeostationaryFactory
from mpop.projector import get_area_def
from pyresample import plot
import numpy as np
from pydecorate import DecoratorAGG
import aggdraw
from trollimage.colormap import rainbow, RainRate
from trollimage.image import Image as trollimage
from PIL import ImageFont, ImageDraw
from pycoast import ContourWriterAGG
import sys
# to get the local time
import ephem
# to get the solar reflectance of the IR_039 channel
from pyspectral.near_infrared_reflectance import Calculator
from pyspectral.solar import (SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
from pyspectral.rsr_reader import RelativeSpectralResponse
from pyorbital.astronomy import sun_zenith_angle
In [2]:
def solartime(observer, sun=ephem.Sun()):
sun.compute(observer)
# sidereal time == ra (right ascension) is the highest point (noon)
hour_angle = observer.sidereal_time() - sun.ra
return ephem.hours(hour_angle + ephem.hours('12:00')).norm # norm for 24h
In [ ]:
In [3]:
# set the date
time_slot=datetime.datetime(2017,6,16,12,0,0,0)
#datetime1=datetime.datetime(2017,6,9,12,0)
In [4]:
# set up to be able to calculate solar time
o = ephem.Observer()
o.date = time_slot # some utc time
o.date
Out[4]:
In [ ]:
In [5]:
load_radar=True
load_sat=True
load_knmi = True
# RADAR
prop_str='RATE'
# SATELLITE
#channel_sat=['VIS006','VIS008','IR_016','IR_039','WV_062','WV_073','IR_087','IR_097','IR_108','IR_120','IR_134','HRV']
channel_sat=['VIS006','VIS008','IR_016','IR_039','WV_062','WV_073','IR_087','IR_097','IR_108','IR_120','IR_134']
#KNMI
# ATTENTION: cloud mask channel (cldmask) missing for now because it doesn't have the attribute units
#& I do not have the rights to try & add them
#channel_knmi = ['azidiff','cth','cot','cph','ctt','cwp','dcld','dcot','dcwp','dndv','dreff','qa','reff','satz','sds','sds_cs','sds_diff','sds_diff_cs','sunz']
# here without sunz because it differs (up to 1.18° in time instance I tried) from the pytroll computed one
#channel_knmi = ['azidiff','cth','cot','cph','ctt','cwp','dcld','dcot','dcwp','dndv','dreff','qa','reff','satz','sds','sds_cs','sds_diff','sds_diff_cs']
# test without cldmask, cph & qa because they give the following mistake:
# TypeError: Fill value 9.96920996839e+36 overflows dtype int16
channel_knmi = ['azidiff','cth','cldmask','cot','cph','ctt','cwp','dcld','dcot','dcwp','dndv','dreff','qa','reff','satz','sds','sds_cs','sds_diff','sds_diff_cs']
#channel_knmi = ['azidiff','cth','cot','ctt','cwp','dcld','dcot','dcwp','dndv','dreff','reff','satz','sds','sds_cs','sds_diff','sds_diff_cs']
In [6]:
if load_radar:
global_radar = GeostationaryFactory.create_scene("odyssey", "", "radar", time_slot)
global_radar.load([prop_str])
print(global_radar)
print("=========================")
if load_sat:
global_sat = GeostationaryFactory.create_scene("meteosat", "09", "seviri", time_slot)
#global_sat = GeostationaryFactory.create_scene("Meteosat-9", "", "seviri", time_slot)
# test to be able to apply get_reflectance -> doesn't work -.- need workaround -> only change name
# right when apply it
#global_sat.load(['IR_108'], reader_level="seviri-level2")
global_sat.load(channel_sat, reader_level="seviri-level2")
print(global_sat)
print("=========================")
if load_knmi:
global_knmi = GeostationaryFactory.create_scene("cpp", "10", "seviri", time_slot)
global_knmi.load(channel_knmi)
print(global_knmi)
print("=========================")
In [ ]:
In [7]:
area='EuropeOdyssey95'
reproject=True
if reproject:
print ("start projection")
# PROJECT data to new area
if load_radar:
data_radar = global_radar.project(area, precompute=True)
#data[prop_str].product_name = global_radar[prop_str].product_name
#data[prop_str].units = global_radar[prop_str].units
global_radar = data_radar
if load_sat:
data_sat = global_sat.project(area, precompute=True)
global_sat = data_sat
from pyorbital.astronomy import sun_zenith_angle
lonlats = global_sat[channel_sat[0]].area.get_lonlats()
#sza = sun_zenith_angle(datetime1, lonlats[0], lonlats[1])
sza = sun_zenith_angle(time_slot, lonlats[0], lonlats[1])
if load_knmi:
data_knmi = global_knmi.project(area, precompute=True)
global_knmi = data_knmi
print('ok')
In [ ]:
In [8]:
# correct the VIS channels for the solar zenith angle
cos_sza = np.cos(np.deg2rad(sza))
global_sat['VIS006']=global_sat['VIS006']/cos_sza
global_sat['VIS008']=global_sat['VIS008']/cos_sza
global_sat['IR_016']=global_sat['IR_016']/cos_sza
In [9]:
#carry out a parallax correction for all channels
global_sat = global_sat.parallax_corr(fill='bilinear', estimate_cth='standard', replace=True)
In [ ]:
In [10]:
# estimate the solar reflectance of the IR_039 chanel (ATTENTION: not trustworthy!)
# If I use the method from the demo program get slightly different values (& since absolute values so
# tiny, this might actually have an impact)
# what should I do with this channel?!
global_sat['IR_039'].info['satname']='Meteosat-9'
global_sat['IR_039'].info['satnumber']=''
refl039=global_sat['IR_039'].get_reflectance(tb11=global_sat['IR_108'].data,sun_zenith=sza,tb13_4=global_sat['IR_134'].data)
# is this one sunzen corrected yet?!?! if not: can just do after/before parallax corr, right?
# so far assume it is... considering I need to provide sza as an input argument ^^
In [ ]:
In [ ]:
In [11]:
area_size = global_sat["IR_108"].data.size
write_netCDF=True
if write_netCDF:
#nc_outfile = '/data/COALITION2/PicturesSatellite/results_BEL/data/logreg/ODYRATE_SEVIRI_' + datetime1.strftime("%Y%m%d_%H%M") + '.nc'
nc_outfile = '/data/COALITION2/PicturesSatellite/results_BEL/data/logreg/ODYRATE_SEVIRI_' + time_slot.strftime("%Y%m%d_%H%M") + '.nc'
print("*** write data to: ncdump ", nc_outfile)
import netCDF4 as nc4
f = nc4.Dataset(nc_outfile,'w', format='NETCDF4')
tempgrp = f.createGroup('collocated odyssey rain rate and SEVIRI reflectances and brightness temperatures')
m = global_radar[prop_str+'-MASK'].data.reshape(area_size)
y = global_radar[prop_str].data.reshape(area_size)
y = y[m==False]
print("dataset length", len(y))
#xx = global_sat['IR_108'].data.reshape(area_size)
#xx = xx[m==False]
#import pdb
#pdb.set_trace()
lon = lonlats[0].reshape(area_size)
lat = lonlats[1].reshape(area_size)
sza = sza.reshape(area_size)
refl039 = refl039.reshape(area_size)
lon = lon[m==False]
lon_rad = np.deg2rad(lon)
lat = lat[m==False]
sza = sza[m==False]
refl039 = refl039[m==False]
cos_sza = np.cos(np.deg2rad(sza))
lstime = np.empty(lon_rad.shape) #local solar time
for i in range(len(lon_rad)):
o.lon = lon_rad[i]
lstime[i]=solartime(o) # = float number that represents an angle in radians and converts to/from a string
#for i in range(100):
# print y.tolist(0)[i], xx[i], lon[i], lat[i], np.cos(np.deg2rad(sza[i]))
#tempgrp.createDimension('idata', area_size)
tempgrp.createDimension('ndata', len(y))
#IR_108 = tempgrp.createVariable('IR_108', 'f4', 'idata')
varids_sat = range(len(channel_sat))
for chn,i in zip(channel_sat,range(len(channel_sat))):
print chn, i
varids_sat[i] = tempgrp.createVariable(chn, 'f4', 'ndata')
varids_knmi =range(len(channel_knmi))
for chn,i in zip(channel_knmi,range(len(channel_knmi))):
print chn, i
varids_knmi[i] = tempgrp.createVariable(chn, 'f4', 'ndata')
ODY_RAINRATE = tempgrp.createVariable('ODY_RAINRATE','f4', 'ndata')
#time = tempgrp.createVariable('Time', 'i4', 'time')
longitude = tempgrp.createVariable('Longitude', 'f4', 'ndata')
latitude = tempgrp.createVariable('Latitude', 'f4', 'ndata')
csza = tempgrp.createVariable('Cosine Solar Zenith Angle','f4', 'ndata')
lst = tempgrp.createVariable('Local Solar Time','f4', 'ndata')
r39 = tempgrp.createVariable('IR_039_refl','f4', 'ndata')
ODY_RAINRATE[:] = y.tolist(0)
longitude[:] = lon
latitude[:] = lat
csza[:] = cos_sza
lst[:] = lstime
r39[:] = refl039
for chn,i in zip(channel_sat,range(len(channel_sat))):
x = global_sat[chn].data.reshape(area_size)
x = x[m==False]
print chn, i
varids_sat[i][:] = x
for chn,i in zip(channel_knmi,range(len(channel_knmi))):
x = global_knmi[chn].data.reshape(area_size)
x = x[m==False]
print chn, i
varids_knmi[i][:] = x
#IR_108[:] = x.tolist(0)
#IR_108[:] = global_sat["IR_108"].data.reshape(area_size) #The "[:]" at the end of the variable instance is necessary
f.close()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: