%Michael Angus 2020 %Quantile mapping, a general script designed for inputing types of data and %choosing the method. Will then do whatever you want based on the method. %First load the forecast and observations you want to use for this case. %Then put the two matricies as variables to include in the loop. %First add paths to forecast 12km and IMD data. addpath 'C:\Users\micha\Documents\Work_HEPPI\regridding_martin' addpath 'C:\Users\micha\Documents\Work_HEPPI\Data\IMD\IMD_conv' load IMD_processed load fr12km_interp_IMD %% %need to rotate the IMERG input in the second dimension. fcst_mat=permute(flip(fr12km_to_IMD,2),[2 3 1]); %Aphrodite is single, and needs to be in double point format. obs_mat=IMD_map_rorder; QM_mat=zeros(size(fcst_mat)); %for memory purposes clear all the bits we don't need. clear IMD_map_rorder fr12km_to_IMD %% choose method. %options: %Gamma distribution 'Gamma' %Generalised extreme value theory 'GEV' %Pasten-Zapata Double Gamma distribution 'DGQM', from Pasten-Zapata et al. %2020 %GEV for 90% up, gamma for 90% down 'GEVG' %Empirical, actual values. 'Emp' meth='DGQM'; %Begin Quantile mapping approach. %for each grid cell, take the points above zero and fit a gamma %distribution. %As a side note, at some point I will need to write some sort of script to %work out what days and ensemble members the forecast data actually is %because you've smashed it all together into one big data set. for i=1:size(obs_mat,1) for j=1:size(obs_mat,2) %for both, take the points in i and j that are above zero and index %the locations. obs_temp=squeeze(obs_mat(i,j,:)); obs_fit_points=find(obs_temp>0); if isempty(obs_fit_points) QM_mat(i,j,:)=NaN; else obs_vals=obs_temp(obs_fit_points); %repeat for forecast fcst_temp=squeeze(fcst_mat(i,j,:)); fcst_fit_points=find(fcst_temp>0); fcst_vals=fcst_temp(fcst_fit_points); %choose distribution here if strcmp(meth,'Gamma') %fit gamma distribution to obs. obs_prms=gamfit(obs_vals); %fit gamma distribution to Forecast fcst_prms=gamfit(fcst_vals); %Fit values for each point based on forecast cdfval=gamcdf(fcst_vals,fcst_prms(1),fcst_prms(2)); %can't cope with exactly 1, make biggest possible. cdfval(cdfval==1)=.9999999999999999; %replace with quantile mapped values. fcst_fitted=gaminv(cdfval,obs_prms(1),obs_prms(2)); elseif strcmp(meth,'GEV') obs_prms=gevfit(obs_vals); fcst_prms=gevfit(fcst_vals); %third parameter required for GEV dist. cdfval=gevcdf(fcst_vals,fcst_prms(1),fcst_prms(2),fcst_prms(3)); fcst_fitted=gevinv(cdfval,obs_prms(1),obs_prms(2),obs_prms(3)); elseif strcmp(meth,'DGQM') %first fit high and low percentiles, at 90% obs_ptile=prctile(obs_vals,90); obs_low=obs_vals(obs_vals=obs_ptile); %repeat for forecast fcst_ptile=prctile(fcst_vals,90); fcst_low=fcst_vals(fcst_vals=fcst_ptile); fcst_fit_low=find(fcst_temp>0 & fcst_temp=fcst_ptile); fcst_fit_points=cat(1,fcst_fit_low,fcst_fit_high); %fit distributions obs_prms_low=gamfit(obs_low); obs_prms_high=gamfit(obs_high); fcst_prms_low=gamfit(fcst_low); fcst_prms_high=gamfit(fcst_high); cdfval_high=gamcdf(fcst_high,fcst_prms_high(1),fcst_prms_high(2)); %can't cope with exactly 1 (inf), make biggest possible. cdfval_high(cdfval_high==1)=.9999999999999999; cdfval_low=gamcdf(fcst_low,fcst_prms_low(1),fcst_prms_low(2)); fcst_fitted_high=gaminv(cdfval_high,obs_prms_high(1),obs_prms_high(2)); fcst_fitted_low=gaminv(cdfval_low,obs_prms_low(1),obs_prms_low(2)); fcst_fitted=cat(1,fcst_fitted_low,fcst_fitted_high); elseif strcmp(meth,'GEVG') %first fit high and low percentiles, at 90% obs_ptile=prctile(obs_vals,90); obs_low=obs_vals(obs_vals=obs_ptile); %repeat for forecast fcst_low=fcst_vals(fcst_vals=obs_ptile); fcst_fit_low=find(fcst_temp>0 & fcst_temp=obs_ptile); fcst_fit_points=cat(1,fcst_fit_low,fcst_fit_high); %fit distributions obs_prms_low=gamfit(obs_low); obs_prms_high=gevfit(obs_high); fcst_prms_low=gamfit(fcst_low); cdfval_low=gamcdf(fcst_low,fcst_prms_low(1),fcst_prms_low(2)); fcst_fitted_low=gaminv(cdfval_low,obs_prms_low(1),obs_prms_low(2)); if isempty(fcst_fit_high) fcst_fitted=fcst_fitted_low; else fcst_prms_high=gevfit(fcst_high); cdfval_high=gevcdf(fcst_high,fcst_prms_high(1),fcst_prms_high(2)); fcst_fitted_high=gevinv(cdfval_high,obs_prms_high(1),obs_prms_high(2)); fcst_fitted=cat(1,fcst_fitted_low,fcst_fitted_high); end elseif strcmp(meth,'Emp') [fcst_ecdf, fcst_ecdf_rain]=ecdf(fcst_vals); [obs_ecdf, obs_ecdf_rain]=ecdf(obs_vals); %for each forecast value, find the threshold it's %closest too from the forecast ecdf. %then from there choose that ecdf value. %Match that to the closest observed ecdf value. %Choose the actual value from that ecdf. fcst_fitted=zeros(size(fcst_vals)); for x=1:length(fcst_vals) [~,y]=min(abs(fcst_vals(x)-fcst_ecdf_rain)); [~,y2]=min(abs(fcst_ecdf(y)-obs_ecdf)); fcst_fitted(x)=obs_ecdf_rain(y2); end else error('method choice not recognized') end %add back to where non zero points are. QM_mat(i,j,fcst_fit_points)=fcst_fitted; end end i end