% IMD rainfall thresholding
% 16-apr-08; AWR

clear
%% load IMD data
% SOURCES .IMD .NCC1-2005 .v2p0 .rf
%   T (Jun-Sep) VALUES
ncload('/Users/andy/data/India/IMD/data2/IMD1degV2_JJAS51-04.nc');
%ncload('/Users/andy/data/India/IMD/data2/IMD1degV1_JJAS51-03.nc');
nlat=length(Y); nlon=length(X); ndat=length(T);

%% mask by selected state

%  SOURCES .IMD .NCC1-2005 .v2p0 .rf
%   T (last) VALUE
%   T removeGRID
%   SOURCES .Features .Political .World .1st_Order .the_geom
%    objectid (769) (841) (863) (902) (915) (916) (925) (937) (953) (962) (969) (983) (995) (999) (1003) (1011) (1016) (1033) (1035) (1036) (1054) (1064) (1094) (1097) (1105) (1130) (1194) (1319) (1352) (1407) (2148) (2457) VALUES
%    rasterize
%   0.5 maskle
%   dataflag
%   0 maskle

ncload('/Users/andy/data/India/India_States_mask.nc')
% state identifiers
% 769	India	Jammu & Kashmir
% 841	India	Himachal Pradesh
% 863	India	Punjab
% 902	India	Uttar Pradesh
% 915	India	Haryana
% 916	India	Chandigarh
% 925	India	Rajasthan
% 937	India	Arunachal Pradesh
% 953	India	Delhi
% 962	India	Sikkim
% 969	India	Assam
% 983	India	Bihar
% 995	India	West Bengal
% 999	India	Nagaland
% 1003	India	Madhya Pradesh
% 1011	India	Meghalaya
% 1016	India	Manipur
% 1033	India	Gujarat
% 1035	India	Tripura
% 1036	India	Mizoram
% 1054	India	Orissa
% 1064	India	Maharashtra
% 1094	India	Daman & Diu
% 1097	India	Dadra and Nagar Haveli
% 1105	India	Andhra Pradesh
% 1130	India	Karnataka
% 1194	India	Goa
% 1319	India	Tamil Nadu
% 1352	India	Kerala
% 1407	India	Pondicherry
% 2148	India	Lakshadweep
% 2457	India	Andaman and Nicobar Islands


% select Orissa (1054)
state=find(objectid==1054);
% select mask
smask=squeeze(dataflag(state,:,:));
indx=find(isfinite(smask));
nsta=length(indx);

%% apply mask
clear rfstate
for t=1:ndat
  rfstate(t,:)=rf(t,indx);
end

%% percentiles of IMD rainfall
p90=prctile(rfstate,90);
p99=prctile(rfstate,99);

% find dates where rainfall exceeds the threshold
clear dates stns
[dates,stns]=find(rfstate>=repmat(p90,ndat,1));
for s=1:nsta
  p90dat{s}=dates(find(stns==s));
end
clear dates stns
[dates,stns]=find(rfland>=repmat(p99,ndat,1));
for s=1:nsta
  p99dat{s}=dates(find(stns==s));
end


  
  




