dheraj281
Posts: 1
Joined: Thu Jun 28, 2018 11:08 pm
Occupation: Student

In ozone profile retrieval, I get many values negative. How can we reduce it?

Postby dheraj281 » Thu Jun 28, 2018 11:12 pm

In ozone profile retrieval, I get many values negative. How can we reduce it? my matlab code for this is:


% This program is to find the ozone profile using INSAT3D data and PCA
% algorithm. Here retrieval has been done for OMPS ozone profiles
% For retrieved ozone profile,
% Pressure_levels = [0.1 0.29 0.69 1.42 2.611 4.407 6.95 10.37 14.81 20.4 27.26 35.51 45.29 56.73 69.97 85.18 102.05 122.04 143.84 167.95 194.36 222.94 253.71 286.6 321.5 358.28 396.81 436.95 478.54 521.46 565.54 610.6 656.43 702.73 749.12 795.09 839.95 882.8 922.46 957.44 985.88 1005.43 1013.25];
% For GOES/INSAT PF=[0.1 0.2 0.5 1.0 1.5 2.0 3.0 4.0 5.0 7.0 10.0 15.0 20.0 25.0 30.0 50.0 60.0 70.0 85.0 100.0 115.0 135.0 150.0 200.0 250.0 300.0 350.0 400.0 430.0 475.0 500.0 570.0 620.0 670.0 700.0 780.0 850.0 920.0 950.0 1000.0];

clc;
clear all;
close all;

path='D:\';
reg_coef_5=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_5_to_0'));
reg_coef_10=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_10_to_0'));
reg_coef_15=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_15_to_0'));
reg_coef_20=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_20_to_0'));
reg_coef_25=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_25_to_0'));
reg_coef_30=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_30_to_0'));
reg_coef_35=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_35_to_0'));
reg_coef_40=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_40_to_0'));
reg_coef_45=load(strcat(path,'PCA\INSAT_validation\reg_insat_deg_45_to_0'));
reg_coef_toz=load(strcat(path,'PCA\INSAT_validation\toz_reg_coeff_40lev_liet_1stalgo_insat3d_sigma'));
reg_coef_PC1=load(strcat(path,'PCA\INSAT_validation\reg_coef_PC1_insat_sigma'));
reg_coef_PC2=load(strcat(path,'PCA\INSAT_validation\reg_coef_PC2_insat_sigma'));
reg_coef_Li_ozoProf=load(strcat(path,'PCA\INSAT_validation\reg_coef_Li_ozo_prof_insat_sigma'));

EoFs=load(strcat(path,'PCA\GOES_validation\EoFs')); % EoFs and ozone mean profile are from seebor dataset, so same casn be used for INSAT and INSAT
Mean_ozone_prof=load(strcat(path,'PCA\GOES_validation\ozone_mean_train'));

sig=[0.0 0.0002 0.0005 0.001 0.0015 0.002 0.003 0.0040 0.005 0.007 0.010 0.015 0.020 0.025 0.030 0.050 0.06 0.07 0.085 0.100 0.115 0.135 0.150 0.200 0.250 0.300 0.350 0.400 0.430 0.475 0.500 0.570 0.620 0.670 0.700 0.780 0.850 0.920 0.950 1.0];
Pressure_levels_INSAT=[0.1 0.2 0.5 1.0 1.5 2.0 3.0 4.0 5.0 7.0 10.0 15.0 20.0 25.0 30.0 50.0 60.0 70.0 85.0 100.0 115.0 135.0 150.0 200.0 250.0 300.0 350.0 400.0 430.0 475.0 500.0 570.0 620.0 670.0 700.0 780.0 850.0 920.0 950.0 1000.0];
PI=Pressure_levels_INSAT;

list=textread(strcat(path,'PCA\INSAT3D_data\2016\Jan2016\list.text'),'%s');
nrow1=length(list);
dir=strcat(path,'PCA\INSAT3D_data\2016\Jan2016\');

list_sbuv=textread(strcat(path,'PCA\OMPS_data\Jan2016\list_omps.text'),'%s');
nrow_sbuv=length(list_sbuv);
dir_ia=strcat(path,'PCA\OMPS_data\Jan2016\');

z=0;
for x=1:nrow1
nam=list{x};
fn=[dir nam];
date=str2num(nam(7:8));
time=str2num(nam(17:20));
hr=str2num(nam(17:18)); minute=str2num(nam(19:20));
minu=hr*60+minute;

sur_pres=hdf5read(fn,'FCST_SURF_PRES'); % hPa
sur_tem=hdf5read(fn,'TSurfPhy');
temp_prof=hdf5read(fn,'TAirPhy');
hum_prof=hdf5read(fn,'H2OMMRPhy'); % in g/kg; fill value -999
toz_insat=hdf5read(fn,'totO3Phy');

zenith=hdf5read(fn,'SAT_ZEN');

lat=hdf5read(fn,'Latitude');
lon=hdf5read(fn,'Longitude');
alat=double(lat)*0.01;
along=double(lon)*0.01;

Si=size(alat);

INSAT_BT=zeros(Si(1),Si(2),18);
INSAT_BT(:,:,1)=hdf5read(fn,'CLR_SKY_BT1');
INSAT_BT(:,:,2)=hdf5read(fn,'CLR_SKY_BT2');
INSAT_BT(:,:,3)=hdf5read(fn,'CLR_SKY_BT3');
INSAT_BT(:,:,4)=hdf5read(fn,'CLR_SKY_BT4');
INSAT_BT(:,:,5)=hdf5read(fn,'CLR_SKY_BT5');
INSAT_BT(:,:,6)=hdf5read(fn,'CLR_SKY_BT6');
INSAT_BT(:,:,7)=hdf5read(fn,'CLR_SKY_BT7');
INSAT_BT(:,:,8)=hdf5read(fn,'CLR_SKY_BT8');
INSAT_BT(:,:,9)=hdf5read(fn,'CLR_SKY_BT9');
INSAT_BT(:,:,10)=hdf5read(fn,'CLR_SKY_BT10');
INSAT_BT(:,:,11)=hdf5read(fn,'CLR_SKY_BT11');
INSAT_BT(:,:,12)=hdf5read(fn,'CLR_SKY_BT12');
INSAT_BT(:,:,13)=hdf5read(fn,'CLR_SKY_BT13');
INSAT_BT(:,:,14)=hdf5read(fn,'CLR_SKY_BT14');
INSAT_BT(:,:,15)=hdf5read(fn,'CLR_SKY_BT15');
INSAT_BT(:,:,16)=hdf5read(fn,'CLR_SKY_BT16');
INSAT_BT(:,:,17)=hdf5read(fn,'CLR_SKY_BT17');
INSAT_BT(:,:,18)=hdf5read(fn,'CLR_SKY_BT18');

for xx=1:nrow_sbuv % some data of date is there in previous day file also; xx=1 corresponds to november 30, xx=2 corresponds to dec 1 and so on
nam_ia=list_sbuv{xx};
dt=str2num(nam_ia(57:58));
if dt==date
fil_nme_ia=[dir_ia nam_ia];
% finding matchup data file w.r.t. time
fid=fopen(fil_nme_ia);
B=fscanf(fid, '%s', 76);
M=fscanf(fid, '%f',inf);
%....................................
% To reshape M in matrix
p=0;
for i=1:size(M,1)
if M(i)==2016
p=p+1;
SBUV(p,1:39)=M(i:i+38,1);
end;
end;

%...............................
day_ia=SBUV(:,2);
date_ia=day_ia; % to have dates in January month
time_sec=SBUV(:,3); lt1=SBUV(:,4); ln1=SBUV(:,5);
for k=1:size(SBUV,1)
minu_ia=time_sec(k)/60.0;
if date==date_ia(k);
if abs(minu_ia-minu)<=60.0 && (lt1(k)<=45.0 && lt1(k)>=0.0) && (ln1(k)>=55.0 && ln1(k)<=100.0); % INSAT sounder data range lat: 5N to 38N and lon: 58E to 100E
%t=abs(minu_ia-minu)
mindist=1.0;
for i=1:Si(1)
for j=1:Si(2)
dist1=sqrt((alat(i,j)-lt1(k))^2 + (along(i,j)-ln1(k))^2);
if dist1<mindist && sur_pres(i,j)>0 && temp_prof(i,j,40)>0 % temp_prof>0 is used because some times there is no retrieval at nearest point
mindist=dist1;
xi=i;
xj=j;
end
end
end
if mindist<1.0
zn=zenith(xi);
if zn<=47.5 && zn>0.0
% to get BT at zero zenith
if zn>2.5 && zn<=7.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_5(:,j);
end;
end;

if zn>7.5 && zn<=12.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_10(:,j);
end;
end;

if zn>12.5 && zn<=17.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_15(:,j);
end;
end;

if zn>17.5 && zn<=22.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_20(:,j);
end;
end;

if zn>22.5 && zn<=27.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_25(:,j);
end;
end;

if zn>27.5 && zn<=32.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_30(:,j);
end;
end;

if zn>32.5 && zn<=37.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_35(:,j);
end;
end;

if zn>37.5 && zn<=42.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_40(:,j);
end;
end;

if zn>42.5 && zn<=47.5
for j=1:18
INSAT_BT_0(j)=[1 INSAT_BT(xi,xj,j)*secd(zn)]*reg_coef_45(:,j);
end;
end;
%..........................

% Interpolation of profiles on sigma levels
TemI=temp_prof(xi,xj,:);
HumI=hum_prof(xi,xj,:);
TemF=zeros(1,40);
HumF=zeros(1,40);
nI=40;
nF=40;

Pa=((sur_pres(xi,xj)-0.1)*sig)+0.1; % 0.1 hPa is pressure at top level
PF=Pa;
ij=1;
for i=1:nF
for j=ij:nI
if sur_pres(xi,xj)>=PI(j)
if (PF(i)-PI(j))==0.0
TemF(1,i)=TemI(1,j);
HumF(1,i)=HumI(1,j);
ij=j;
break;
end;
if (PF(i)-PI(j))<0.0
TemF(1,i)=TemI(1,j-1)+(log(PF(i)/PI(j-1))*((TemI(1,j)-TemI(1,j-1))/(log(PI(j)/PI(j-1)))));
HumF(1,i)=HumI(1,j-1)+(log(PF(i)/PI(j-1))*((HumI(1,j)-HumI(1,j-1))/(log(PI(j)/PI(j-1)))));
ij=j;
break;
end;
end;
if sur_pres(xi,xj)<PI(j)
if (PF(i)-PI(j))<0.0
TemF(1,i)=TemI(1,j-2)+(log(PF(i)/PI(j-2))*((TemI(1,j-1)-TemI(1,j-2))/(log(PI(j-1)/PI(j-2)))));
HumF(1,i)=HumI(1,j-2)+(log(PF(i)/PI(j-2))*((HumI(1,j-1)-HumI(1,j-2))/(log(PI(j-1)/PI(j-2)))));
ij=j;
break;
end;
end;
end;
end;
if (PF(nF)-PI(nI))>0.0
TemF(1,nF)=TemI(1,nI-1)+(log(PF(nF)/PI(nI-1))*((TemI(1,nI)-TemI(1,nI-1))/(log(PI(nI)/PI(nI-1)))));
HumF(1,nF)=HumI(1,nI-1)+(log(PF(nF)/PI(nI-1))*((HumI(1,nI)-HumI(1,nI-1))/(log(PI(nI)/PI(nI-1)))));
end;
%.................................

% total ozone calculation
%X=[1 INSAT_BT_0(9) (INSAT_BT_0(9).^2)/250 sur_tem(xi,xj) temp_prof(xi,xj,:) log(hum_prof(xi,xj,23:40))]; % hum is not multiplied by 1000 as it is already in g/kg
X=[1 INSAT_BT_0(9) (INSAT_BT_0(9).^2)/250 sur_tem(xi,xj) TemF log(HumF(1,23:40))]; % hum is not multiplied by 1000 as it is already in g/kg
TOZ=X*reg_coef_toz; % reg_coef_toz is 62 *1

% Finding of PCs
X_PC1=[INSAT_BT_0(3:4) INSAT_BT_0(9) INSAT_BT_0(13:15) 1];
PC1_ret=X_PC1*reg_coef_PC1; % reg_coef_PC1 is 7*1
X_PC2=[TOZ 1];
PC2_ret=X_PC2*reg_coef_PC2;

% To retrieve ozone profile using retrieved PC1_ret and PC2_ret
O3_prof_ret = Mean_ozone_prof + PC1_ret*(EoFs(:,1)') + PC2_ret*(EoFs(:,2)'); % retrieved ozone profile is from top to bottom
%..........................................................................

% Retrieval by Li's method
latitude=cos(alat(xi,xj)*3.14/180)*10.0+200.0;
mon=12; % because this is for december
month=cos((mon-6)*3.14/12)*10.0+200.0;
Z=[INSAT_BT_0(1:15) (INSAT_BT_0(1:15).^2)/250 sur_pres(xi,xj) latitude month 1];
log_Ozo_prof_Li=Z*reg_coef_Li_ozoProf;
Ozo_prof_Li=exp(log_Ozo_prof_Li);

z=z+1;
%A(z,1)=yr;
A(z,2)=minu_ia-minu; A(z,3)=date; A(z,4)=minu; A(z,5)=minu_ia;
A(z,6)=alat(xi,xj); A(z,7)=along(xi,xj); A(z,8)=lt1(k); A(z,9)=ln1(k); A(z,10)=mindist;
A(z,11)=sur_pres(xi,xj);A(z,12)=sur_tem(xi,xj);
%A(z,13)=ind_land_water(xi,xj);A(z,14)=sur_elev(xi,xj);
A(z,15)=zn;
A(z,16:55)=O3_prof_ret; A(z,56:94)=SBUV(k,:);
A(z,95:134)=Ozo_prof_Li; A(z,135)=x; A(z,136)=xx;
%save(strcat(path,'PCA\INSAT_validation\INSAT_SBUV'),'A','-ASCII');
end
end
end
end
end
end;
end
clearvars -EXCEPT reg_coef_5 reg_coef_10 reg_coef_15 reg_coef_20 reg_coef_25 reg_coef_30 reg_coef_35 reg_coef_40 reg_coef_45 reg_coef_toz reg_coef_PC1 reg_coef_PC2 reg_coef_Li_ozoProf EoFs Mean_ozone_prof sig Pressure_levels_INSAT PI list nrow1 dir list_sbuv nrow_sbuv dir_ia A z
fclose('all');
end
save('D:\PCA\INSAT_validation\INSAT_OMPS_Dheeraj_60min','A','-ASCII');

AmyCowen
Site Admin
Posts: 104
Joined: Mon Aug 22, 2016 4:39 pm
Occupation: Administrator

Re: In ozone profile retrieval, I get many values negative. How can we reduce it?

Postby AmyCowen » Thu Jul 19, 2018 5:02 pm

Hi - This forum is for students using the Science Buddies Raspberry Pi Projects Kit.

Science Buddies


Return to “Raspberry Pi Projects Kit”