diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f9cbc9e --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +data/ +prelres/ \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..a8fc851 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2016 + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/OT_d_best_20050311.mat b/OT_d_best_20050311.mat new file mode 100644 index 0000000..4d361ea Binary files /dev/null and b/OT_d_best_20050311.mat differ diff --git a/OT_d_best_20050314.mat b/OT_d_best_20050314.mat new file mode 100644 index 0000000..5147fef Binary files /dev/null and b/OT_d_best_20050314.mat differ diff --git a/OT_models_best.mat b/OT_models_best.mat new file mode 100644 index 0000000..0f18c8d Binary files /dev/null and b/OT_models_best.mat differ diff --git a/README.md b/README.md new file mode 100644 index 0000000..d9e8caf --- /dev/null +++ b/README.md @@ -0,0 +1,8 @@ +# VIEWS +The Lake Victoria Intense storm Early Warning System (VIEWS) + +Version 1.0 of the package, termed Lake Victoria Intense storm Early Warning System (VIEWS), is available at https://github.com/wthiery/VIEWS and is released under the MIT licence. + +At this stage the prediction system needs to be considered as a prototype; more research as well as input from the user community is needed to improve its skill. + +When applied, the package runs permanently, and is activated at forecast lead time. The software first downloads the OT images corresponding to the daytime hours from the NASA Langley Research Centre data repository. It subsequently computes the predictor value OT_{day} for each country and the whole lake by performing the appropriate spatial and temporal selection (see equation 2 and table 1 in Thiery et al., 2017 ERL). The OT_day values then serve as input for the respective logistic regressions (see equation 1 and table 1 in Thiery et al., 2017 ERL), yielding the probability for an extreme event. Depending on the threshold probability defined by the user, the software will indicate whether or not a warning is to be issued for a specific lake sector or the whole lake. diff --git a/main.m b/main.m new file mode 100644 index 0000000..74bba6d --- /dev/null +++ b/main.m @@ -0,0 +1,137 @@ + +% -------------------------------------------------------------------- +% The Lake Victoria Intense storm Early Warning System (VIEWS) +% +% When using this model please cite this reference: +% +% Thiery, W., Gudmundsson, L., Bedka, K., Semazzi, F.H.M., Lhermitte, +% S., Willems, P., van Lipzig, N.P.M. and Seneviratne, S.I. Early +% warnings of hazardous thunderstorms on Lake Victoria, Env. Res. +% Lett., in review. +% -------------------------------------------------------------------- + + +% to do: +% - hours_day differs between whole lake, KEN, TAN and UGA !!! +% - 'clock' should be UTC not local time !!! + + + +% -------------------------------------------------------------------- +% main script to perform intense storm prediction +% tested on MATLAB 7.12.0 (R2011a) +% -------------------------------------------------------------------- + + +% start clock +tic + + +% clean up +clc; +clear; +close all; + + +% flags +flags.user_OTdata = 2; % 0: use historical data + % 1: request input from user about OT data + % 2: operational mode with data from NASA server +flags.user_config = 0; % 0: use optimised model configuration from Thiery et al., 2017 ERL + % 1: request input from user about model configuration +flags.plot = 0; % 0: do not plot + % 1: plot + + + +% -------------------------------------------------------------------- +% initialisation +% -------------------------------------------------------------------- + + +% define hours considered as "day" and "night" +hours_night = [22:24 1:9]; % final estimate based on OT diurnal cycle +hours_day = 10:15; + + +% define percentile above which events are considered 'severe' +perc_severe = 99; + + +% initialise model parameters +res_reg = 0.20; % predefined regular grid + + +% name of the best logistic regression model (obtained by optimisation, see paper section 4) +model_best = 'OT_models_best.mat'; + + +% name of the best logistic regression model (obtained by optimisation, see paper section 4) +% OT_data_test = 'OT_d_best_20050311.mat'; % example data for testing the model - calm day (11/03/2005) +OT_data_test = 'OT_d_best_20050314.mat'; % example data for testing the model - extreme day (14/03/2005) + + +% Define remote path where OT data can be downloaded +OT_rpath = 'https://clouds.larc.nasa.gov/prod/exp/lake_victoria'; + + +% add data path to work path +addpath('data') + + +% -------------------------------------------------------------------- +% load data +% -------------------------------------------------------------------- + + +% load best logistic regression model +if flags.user_config == 0 + load(model_best) +elseif flags.user_config == 1 + model_user = input('Please name the file containing the logistic regression model (e.g. LOGR_OT_best.mat)\n\n','s'); + load(model_user) +end + + +% load Overshooting Top (OT) data +if flags.user_OTdata == 0 % use historical data + + disp(sprintf(['Loading test data (' OT_data_test ')\n'])) %#ok<*DSPS> + load(OT_data_test) + +elseif flags.user_OTdata == 1 % request input from user about OT data + + OT_data_user = input('Please name the file containing the afternoon OT map (e.g. OT_d_best_20050311.mat)\n\n','s'); + load(OT_data_user) + +elseif flags.user_OTdata == 2 % operational mode with data from NASA server + + [OT_d, OT_d_regridded_daysum] = mf_get_OT_today(OT_rpath, LOGR_OT_best.hours_day, LOGR_OT_best.lat, LOGR_OT_best.lon); + +end + + + +% -------------------------------------------------------------------- +% Perform prediction +% -------------------------------------------------------------------- + + +% whole lake +mf_VIEWS('Whole lake', OT_d, LOGR_OT_best); + + +% Sector Uganda +mf_VIEWS('Uganda' , OT_d, LOGR_OT_best_Uga); + + +% Sector Kenya +mf_VIEWS('Kenya' , OT_d, LOGR_OT_best_Ken); + + +% Sector Tanzania +mf_VIEWS('Tanzania' , OT_d, LOGR_OT_best_Tan); + + +% stop clock +toc diff --git a/mf_VIEWS.m b/mf_VIEWS.m new file mode 100644 index 0000000..c1ceb03 --- /dev/null +++ b/mf_VIEWS.m @@ -0,0 +1,37 @@ + + +% -------------------------------------------------------------------- +% function to perform logistic regression and compute ROC curve +% -------------------------------------------------------------------- + + +function [warning] = mf_VIEWS(sector, OT_d, LOGR_OT) + + +% -------------------------------------------------------------------- +% manipulations +% -------------------------------------------------------------------- + + +% get daytime daytime total OTs over highly correlated regions +OT_d_gp = nansum(nansum(OT_d(LOGR_OT.iscorr),1),2); + + +% get probability for an extreme night given the input logistic regression model +p_all = glmval(LOGR_OT.b_all, OT_d_gp, 'logit'); + + +% decide whether warning nees to be issued +% assume optimal point (threshold probability for the difference between +% hit rate and false alarm rate is maximum) +if p_all < LOGR_OT.T_opt % do not issue warning + warning = 0; + disp(sprintf(['NO WARNING (' sector ')\n'])) %#ok<*DSPS> +elseif p_all >= LOGR_OT.T_opt % issue warning + warning = 1; + disp(sprintf(['WARNING (' sector '):\n High probability for extreme nighttime thunderstorm activity on Lake Victoria\n'])) +end + + + +end diff --git a/mf_get_OT_today.m b/mf_get_OT_today.m new file mode 100644 index 0000000..9550309 --- /dev/null +++ b/mf_get_OT_today.m @@ -0,0 +1,197 @@ + + +% -------------------------------------------------------------------- +% function to perform logistic regression and compute ROC curve +% -------------------------------------------------------------------- + + +function [OT_d_regridded, OT_d_regridded_daysum] = mf_get_OT_today(OT_rpath, hours_day, lat, lon) + + + +% flags +flags.archive = 1; % 0: Do not archive OT data (i.e. delete files) + % 1: Archive OT data (i.e. move files to archive directory) + + +% -------------------------------------------------------------------- +% initialisation +% -------------------------------------------------------------------- + + +% initialise basename of the file +basename = 'NASA_LARC_SEVIRI_OTDETECTION_'; + + +% initialise grid size from original msg images +nlon_msg = 560; +nlat_msg = 672; + + + +% -------------------------------------------------------------------- +% manipulations +% -------------------------------------------------------------------- + + +% get time info +now.date_vec = clock; +now.year = now.date_vec(1); +now.month = now.date_vec(2); +now.day = now.date_vec(3); +now.hour = now.date_vec(4); +now.min = now.date_vec(5); +now.sec = now.date_vec(6); +now.dayofyr = nansum(eomday(now.year, 1:now.month-1)) + now.day; +disp(sprintf(['Time is ' num2str(now.hour, '%02d') ':' num2str(now.min, '%02d') 'h on ' num2str(now.day) '/' num2str(now.month) '/' num2str(now.year) '\n'])) %#ok<*DSPS> + + +% % % % % TESTING +% % % for ind = 97:110; +% % % now.dayofyr = ind; +now.dayofyr = 036; +now.hour = 18; +% % % % % TESTING + + +% check that prediction is made during the richt time of the day +if now.hour < hours_day(end) + error('myApp:timeofday', ['warnings can only be delivered between ' num2str(hours_day(end)) 'h and 23.59h UTC']) +end + + +% construct output file directory, name and url +filedir = [num2str(now.year, '%04d') num2str(now.dayofyr, '%03d')]; + + +% create archive directory +if flags.archive == 1 + archivedir = ['archive/' filedir]; + mkdir(archivedir); +end + + +% print status message to screen +disp(sprintf(['downloading OT data from ' OT_rpath])) + + +% initialise variable +OT_d = 0; + + +% loop over files +for i=1:length(hours_day); + for j=[00 15 30 45] + + + % construct output file name and url + filename = [basename num2str(now.year, '%04d') num2str(now.dayofyr, '%03d') '.' num2str(hours_day(i), '%02d') num2str(j, '%02d') '.nc']; + url = [OT_rpath '/' filedir '/' filename]; + + + % download OT data from remote path only if file is there + [status, ~] = system(['wget -P data/ -N ' url]); + + + % if dowload was succesful: load and process data + if status == 0 % download succesful + + + % print status message to screen + disp(sprintf([filename ' succesfully downloaded'])) + + + % fix error with old matlab versions, this should not be used anymore in newer matlab versions + %mf_fix_netcdf4_dimid(filename); + + + % load raw variables + lat_msg = ncread(filename, 'latitude' ) ; + lon_msg = ncread(filename, 'longitude' ) ; + ir_brightness_temperature = rot90(ncread(filename, 'ir_brightness_temperature' )); + try + tropopause_temperature = rot90(ncread(filename, 'tropopause_temperature' )); + ot_anvilmean_brightn_temp_diff = rot90(ncread(filename, 'ot_anvilmean_brightness_temperature_difference')); + catch errmessage %#ok<*NASGU> + disp('failed to load all required variables; assume no OTs were detected on this image') + tropopause_temperature = NaN(size(ir_brightness_temperature)); + ot_anvilmean_brightn_temp_diff = NaN(size(ir_brightness_temperature)); + end + + + % detect OT from raw variables, by looking for pixels meeting each of the following 3 conditions: + cond1 = ir_brightness_temperature < 217.5; % colder than 217.5 K (ir_brightness_temperature in NetCDF file) + cond2 = ir_brightness_temperature - tropopause_temperature <= 5; % ir_brightness_temperature up to 5 K warmer than the tropopause temperature (tropopause_temperature in NetCDF file) + cond3 = ot_anvilmean_brightn_temp_diff >= 6; % ir_brightness_temperature 6 K or more colder than the anvil (ot_anvilmean_brightness_temperature_difference in NetCDF file) + OT_15min = cond1 .* cond2 .* cond3; + + + % if dowload was NOT succesful: assume no OTs were detected + else % dowload unsuccesful + + + disp(sprintf(['failed to download ' filename '; assume no OTs were detected'])) + OT_15min = zeros(nlat_msg, nlon_msg); + + + end + + + % aggregate all OT + OT_d = OT_d + OT_15min; + + +% % archive or delete OT data +% if flags.archive == 0 +% delete(filename); % delete file +% elseif flags.archive == 1 +% movefile(filename,archivedir,'f') % copy file +% end + + + end + +end + + +% debugging +% figure;imagesc(OT_d);colorbar +% debugging + + +% print white line to the screen +disp(sprintf('\n')) + + +% brute-force regridding: +[nlat, nlon] = size(lat); +OT_d_regridded = zeros(nlat, nlon); +res_lat = abs( lat(2,1) - lat(1,1) ); +res_lon = abs( lon(1,2) - lon(1,1) ); + + +% loop over pixels in the new coarse grid +for i=1:nlat + for j=1:nlon + + % get the corresponding pixels in the fine grid + rows_msg = find(lat_msg >= lat(i,1)-res_lat/2 & lat_msg < lat(i,1)+res_lat/2); + cols_msg = find(lon_msg >= lon(1,j)-res_lon/2 & lon_msg < lon(1,j)+res_lon/2); + + % sum them, and assign the coarse pixel this value + OT_d_regridded(i,j) = nansum(nansum(OT_d(rows_msg, cols_msg))); + + end +end + + +% % % % debugging +% % % figure;imagesc(OT_d_regridded);colorbar; caxis([0 200]); +% % % OT_d_regridded_daysum(ind) = nansum(nansum(OT_d_regridded)); +OT_d_regridded_daysum = NaN; +% % % % debugging + +% % % end + + +end diff --git a/wget.exe b/wget.exe new file mode 100644 index 0000000..2d2898d Binary files /dev/null and b/wget.exe differ