In [1]:
%matplotlib inline
import sys
import time
from datetime import datetime
sys.path.append('..')
from netCDF4 import Dataset
import pylab as pl
import numpy as np
import scipy.interpolate as interp
%matplotlib inline
def cosd(inn):
return np.cos(inn*np.pi/180.)
def sind(inn):
return np.sin(inn*np.pi/180.)
def acosd(inn):
return np.arccos(inn)*180./np.pi
def asind(inn):
return np.arcsin(inn)*180./np.pi
def azi2azid(sa,va):
return acosd(cosd(sa)*cosd(va)+sind(sa)*sind(va))
def height2press(hh):
return 1013.*(1.-(hh*0.0065/288.15))**5.2555
def histeq(im,nbr_bins=256):
imhist,bins = np.histogram(im.flatten(),nbr_bins,normed=True)
cdf = imhist.cumsum() #cumulative distribution function
cdf = 255 * cdf / cdf[-1] #normalize
im2 = np.interp(im.flatten(),bins[:-1],cdf)
return im2.reshape(im.shape)
def get_altitude(ff,st=(1,1)):
with Dataset('%s/geo_coordinates.nc'%ff) as ds:
alt=ds.variables['altitude'][::st[0],::st[1]]*1.
return alt
def interpolate_tie(inn,al,ac):
npa=np.arange
npl=np.linspace
sh_in=inn.shape
sh_ou=((sh_in[0]-1)*al+1,(sh_in[1]-1)*ac+1)
ou=interp.RectBivariateSpline(
npl(0.,sh_ou[0],sh_in[0])
, npl(0.,sh_ou[1],sh_in[1])
, inn, kx=1, ky=1)(npa(sh_ou[0])
, npa(sh_ou[1]))
return ou
def get_pressure(ff,st=(1,1)):
with Dataset('%s/tie_meteo.nc'%ff) as ds:
lat=ds.variables['sea_level_pressure'][:]*1.
ac=ds.ac_subsampling_factor
al=ds.al_subsampling_factor
return interpolate_tie(lat,al,ac)[::st[0],::st[1]]
def get_latitude(ff,st=(1,1)):
## assumes FR dataset
with Dataset('%s/tie_geo_coordinates.nc'%ff) as ds:
lat=ds.variables['latitude'][:]*1.
ac=ds.ac_subsampling_factor
al=ds.al_subsampling_factor
return interpolate_tie(lat,al,ac)[::st[0],::st[1]]
def get_longitude(ff,st=(1,1)):
## assumes FR dataset
with Dataset('%s/tie_geo_coordinates.nc'%ff) as ds:
lon=ds.variables['longitude'][:]*1.
ac=ds.ac_subsampling_factor
al=ds.al_subsampling_factor
re=np.sin(lon*np.pi/180.)
im=np.cos(lon*np.pi/180.)
re_out=interpolate_tie(re,al,ac)[::st[0],::st[1]]
im_out=interpolate_tie(im,al,ac)[::st[0],::st[1]]
return np.arctan2(re_out,im_out)*180./np.pi
def get_geometry(ff,st=(1,1)):
out={}
with Dataset('%s/tie_geometries.nc'%ff) as ds:
ac=ds.ac_subsampling_factor
al=ds.al_subsampling_factor
out={k: ds.variables[k][:]*1. for k in ds.variables}
for k in out:
out[k]=interpolate_tie(out[k],al,ac)[::st[0],::st[1]]
out['ADA']=azi2azid(out['SAA'],out['OAA'])
return out
def get_cwl(ff,bn,st=(1,1)):
with Dataset('%s/instrument_data.nc'%ff) as ds:
l0=(ds.variables['lambda0'][bn-1,:]*1.).squeeze()
di=ds.variables['detector_index'][:]*1
return l0[di][::st[0],::st[1]]
def get_solar_flux(ff,bn,st=(1,1)):
with Dataset('%s/instrument_data.nc'%ff) as ds:
sf=(ds.variables['solar_flux'][bn-1,:]*1.).squeeze()
di=ds.variables['detector_index'][:]*1
return sf[di][::st[0],::st[1]]
def get_albedo(lat,lon,doy):
# this is meris white sky surface albedo of band 10
# from albedomap year 2005
# better data e.g. from CAWA aerosol is welcome :)
lat_idx=np.round((90. -lat)*20).astype(np.int).clip(0,3599)
lon_idx=np.round((180.+lon)*20).astype(np.int).clip(0,7199)
with Dataset('../luts/ws_alb_10_2005.nc') as nc:
#get closest day of year
doy_idx=np.abs(nc.variables['time'][:]-doy).argmin()
alb = nc.variables['albedo'][doy_idx,:,:]
#get closest albedo
#nearest neighbour
return alb[lat_idx,lon_idx].clip(0,1.)
def get_doy(ff):
with Dataset('%s/time_coordinates.nc'%ff) as ds:
doy=datetime.strptime(ds.start_time,"%Y-%m-%dT%H:%M:%S.%fZ").timetuple().tm_yday
return doy
def get_band(ff,bn,st=(1,1)):
variable='Oa%02i_radiance'%bn
with Dataset('%s/%s.nc'%(ff,variable)) as ds:
out=ds.variables[variable][:]*1.
return out[::st[0],::st[1]]
olci='S3A_OL_1_EFR____20160531T093141_20160531T093441_20160531T113517_0179_004_364_2160_MAR_O_NR_001.SEN3/'
#pl.imshow(get_band(olci,21))
#pl.imshow(get_cwl(olci,21))
#pl.imshow(get_longitude(olci,))
#pl.imshow(get_geometry(olci)['ADA'])
#pl.imshow(get_altitude(olci).clip(-2500,2500))
#pl.imshow(get_albedo(get_latitude(olci,(20,20)),get_longitude(olci,(20,20))))
#pl.imshow(get_longitude(olci))
#get_albedo(get_latitude(olci,(20,20)),get_longitude(olci,(20,20)))
#pl.imshow(np.fromfile('ws_albedo_10.img',dtype=np.float32).reshape(3600,7200).byteswap())
#pl.colorbar()
get_doy(olci)
Out[1]:
In [2]:
stride=(24,24)
rad={k:get_band(olci,k,stride)/get_solar_flux(olci,k,stride) for k in [2,4,6,12,13,14,15,]}
geo=get_geometry(olci,stride)
lam=get_cwl(olci,13,stride)
lon=get_longitude(olci,stride)
lat=get_latitude(olci,stride)
doy=get_doy(olci)
alb=get_albedo(lat,lon,doy)
prs=height2press(get_altitude(olci,stride))
In [3]:
for xxx in (lam,lon,lat,alb,prs):
print xxx.shape,
print
for xxx in rad:
print rad[xxx].shape,
print
for xxx in geo:
print geo[xxx].shape,
In [4]:
rgb=np.dstack((rad[6],rad[4],rad[2]))
rgb.shape
pl.imshow(rgb/0.3)
pl.imshow(alb)
Out[4]:
In [5]:
import cloud_core
reload(cloud_core)
cc=cloud_core.cloud_core('../luts/cloud_core_olci.nc4')#,used_ab='15')
nnn=rad[12].size
sha=rad[12].shape
oe_ctp=np.zeros(sha)*np.nan
oe_cot=np.zeros(sha)*np.nan
oe_cnv=np.zeros(sha).astype(np.bool)
oe_avk=np.zeros(sha)+1.
oe_cst=np.zeros(sha)+1.
print time.asctime()
a=time.time()
print nnn
jl=0
jw=0l
for i in range(0,nnn):
inp={'suz':geo['SZA'].flat[i]
,'vie':geo['OZA'].flat[i]
,'azi':180. - geo['ADA'].flat[i] # momo definition
,'prs':prs.flat[i] # not yet used!!!!
,'dwl':lam.flat[i]-cc.cha['13']['cwvl'] # first abs-band
,'alb':alb.flat[i] # last win-band
,'rtoa':{'12':rad[12].flat[i]
,'13':rad[13].flat[i]
,'14':rad[14].flat[i]
,'15':rad[15].flat[i]
}
}
erg=cc.estimator(inp)
oe_ctp.flat[i]=erg['ctp']
oe_cot.flat[i]=erg['cot']
oe_cnv.flat[i]=erg['res'].conv
oe_avk.flat[i]=erg['res'].a.trace()
oe_cst.flat[i]=erg['res'].sr[0,0]
jl+=1
if i % 10000 ==0: print i,
print ''
print time.asctime()
In [6]:
pl.figure(figsize=(32,8))
pl.subplot(1,4,2)
pl.imshow(np.where(oe_cot>0.7,oe_ctp,np.nan),aspect='auto')
#pl.imshow(oe_ctp,aspect='auto')
pl.colorbar()
pl.subplot(1,4,3)
pl.imshow(np.where(oe_cot > 0.8,oe_cot,np.nan),aspect='auto')
#pl.imshow(oe_cot,aspect='auto')
#pl.colorbar()
pl.subplot(1,4,1)
rrr=np.dstack((rad[6],rad[4],rad[2]))
rrr.shape
#pl.imshow(rgb[4000//stride[0]:6000//stride[0],:,:]/rgb[:,:,2].max())
pl.imshow(rrr/0.3,aspect='auto')
#pl.subplot(1,4,4)
#pl.imshow(np.where(oe_cot > 0.8,rad['lam11'],np.nan),aspect='auto')
#pl.colorbar()
Out[6]:
In [7]:
pl.figure(figsize=(32,8))
pl.subplot(1,4,2)
pl.imshow(np.where(oe_cot>0.7,oe_ctp,np.nan),aspect='auto')
#pl.imshow(oe_ctp,aspect='auto')
pl.colorbar()
pl.subplot(1,4,3)
pl.imshow(np.where(oe_cot > 0.8,oe_cot,np.nan),aspect='auto')
#pl.imshow(oe_cot,aspect='auto')
#pl.colorbar()
pl.subplot(1,4,1)
rrr=np.dstack((rad[6],rad[4],rad[2]))
rrr.shape
#pl.imshow(rgb[4000//stride[0]:6000//stride[0],:,:]/rgb[:,:,2].max())
pl.imshow(rrr/0.3,aspect='auto')
#pl.subplot(1,4,4)
#pl.imshow(np.where(oe_cot > 0.8,rad['lam11'],np.nan),aspect='auto')
#pl.colorbar()
Out[7]:
In [9]:
#pl.figure(figsize=(32*1.5,24*1.5))
pl.imshow(np.where(oe_cot>0.7,oe_ctp,np.nan))#,aspect='auto')
#pl.imshow(oe_ctp ,aspect='auto')
#pl.colorbar()
pl.axis('off')
#pl.savefig('/home/rene/Downloads/first_olci_ctp_noaxes.png',bbox_inches='tight', pad_inches=0)
Out[9]:
In [131]:
data={
'lon' :{'data':lon.astype(np.float32) ,
'attributes': {
'long_name':'longitude',
'standard_name': 'longitude',
'units':'degrees_east'
}
},
'lat' :{'data':lat.astype(np.float32),
'attributes': {
'long_name':'latitude',
'standard_name': 'latitude',
'units':'degrees_north'
}
},
'ctp':{'data':oe_ctp.astype(np.float32),
'attributes': {
'long_name':'cloud top pressure',
'units':'hPa',
}
},
'cot':{'data':oe_cot.astype(np.float32),
'attributes': {
'long_name':'loh_10 of cloud optical thickness @750 nm',
'units':'1',
}
}
}
from dict2beamncdf import dict2beamncdf as d2n
ret=d2n(data,'demo_olci_erg.nc')
In [11]: