Calculation of Seasonal Extreme Precipitation Indices
Seasonal Extreme Precipitation Indices
Author: Zhi Li
Email:[email protected]
Author: Xianjun Xiang
Email: [email protected]
1
Initialization Initialization
%% Initialization
clc; % Clear command window
clear; % Clear workspace variables
close all % Close all figure windows
% Define time range
StartYear = 1961; % Starting year of data
EndYear = 2022; % Ending year of data
Year = StartYear : EndYear; % Generate year sequence array2
File Import and Variable ReadingFile import and variable reading
%% File import and variable reading
ID = ['E:\北大兼职\臭氧高温复合\extreme_pre\', 'CHM_PRE_0.1dg_19612022.nc'];
lat = ncread(ID, 'latitude'); % Latitude
lon = ncread(ID, 'longitude'); % Longitude
year = Year - StartYear + 1; % Create year time index
% Extract coordinates of target area
lat = lat(16:86);
lon = lon(366:466);
% Set NetCDF file reading parameters
stlo = 366; stla = 16; stpre = 1; start = [stlo, stla, stpre]; % Starting indices for each dimension
locount = 101; lacount = 71; ticount = 22645;count = [locount, lacount, ticount]; % Number of elements to read
stride = [1, 1, 1]; % Stride for reading data
pre_spei_stride1 = ncread(ID, 'pre', start, count, stride);
pre = double(pre_spei_stride1); % Convert precipitation data to double3
Region Masking Region Masking
%% Region Masking
% Import SHP file of Guangdong province
cdL_NMG=shaperead("E:\北大兼职\臭氧高温复合\extreme_pre\广东省\广东省.shp");
cx_NMG=[cdL_NMG(:).X]; % Longitude
cy_NMG=[cdL_NMG(:).Y]; % Latitude
% Create region mask
for i = 1:101
for j = 1:71 % Verify if point is within Guangdong province
in = inpolygon(lon(i),lat(j),cx_NMG,cy_NMG);
if in == 0 % Outside the boundary
pre(i,j,:) = NaN; % Set to NaN
end
end
end
% Get days in a year
days_in_year = @(y) 365 + (mod(y, 4) == 0 & (mod(y, 100) ~= 0 | mod(y, 400) == 0));
% Initialize extreme precipitation indices
PRCPTOT_seasonal = zeros(locount, lacount, 4, length(Year));
R10_seasonal = zeros(locount, lacount, 4, length(Year)); % R10: Days with precipitation ≥10mm
R20_seasonal = zeros(locount, lacount, 4, length(Year)); % R20: Days with precipitation ≥20mm
% Define seasons by months
seasons = {3:5, 6:8, 9:11, [12, 1, 2]};
% Introduction:
% In natural science research, the seasonal order should follow "winter, spring, summer, autumn." The winter season spans December of the previous year and January–February of the current year.
% For example: Winter 2018 includes December 2017 + January & February 2018
% Attention: For coding simplicity, December 2018 and January–February 2019 are treated as the winter of 2018 in this analysis.4
Main Calculation LoopMain Calculation Loop
% Loop processes each year's data
for y = 1:length(Year) % Calculate total days before current year
days_before = sum(arrayfun(days_in_year, Year(1:y-1))); % Get number of days in current year
num_days = days_in_year(Year(y)); % Get precipitation data for current year
pre_year = pre(:, :, days_before + 1 : days_before + num_days); % Seasonal loop for s = 1:4 % Get months for current season months = seasons{s}; % Initialize seasonal precipitation data pre_season = []; % Special handling for winter season if s == 4 for m = months if m == 12 pre_season = cat(3, pre_season, pre_year(:, :, sum(eomday(Year(y), 1:11)) + 1:sum(eomday(Year(y), 1:12)))); elseif m == 1 || m == 2 % Handle January and February of the next year, check if there is a next year if y < length(Year) days_next_year = days_in_year(Year(y+1)); days_before_next_year = sum(arrayfun(days_in_year, Year(1:y))); pre_next_year = pre(:, :, days_before_next_year + 1 : days_before_next_year + days_next_year); start_day = sum(eomday(Year(y+1), 1:m-1)) + 1; end_day = start_day + eomday(Year(y+1), m) - 1; pre_season = cat(3, pre_season, pre_next_year(:, :, start_day:end_day)); end end end else for m = months start_day = sum(eomday(Year(y), 1:m-1)) + 1; end_day = start_day + eomday(Year(y), m) - 1; pre_season = cat(3, pre_season, pre_year(:, :, start_day:end_day)); end end % Calculate extreme precipitation indices for current season for i = 1:locount for j = 1:lacount PRCPTOT_seasonal(i, j, s, y) = sum(pre_season(i, j, pre_season(i, j, :) > 1), 'all'); R10_seasonal(i, j, s, y) = sum(pre_season(i, j, :) >= 10, 'all'); R20_seasonal(i, j, s, y) = sum(pre_season(i, j, :) >= 20, 'all'); end end end
end5
Masking Data Outside Boundary
Masking Data Outside Boundary
% Masking Data Outside Boundary
for i = 1:101
for j = 1:71
for s = 1:4 % Check if point is within Guangdong in = inpolygon(lon(i), lat(j), cx_NMG, cy_NMG); if in == 0 % Outside the boundary PRCPTOT_seasonal(i, j, s, :) = NaN; R10_seasonal(i, j, s, :) = NaN; R20_seasonal(i, j, s, :) = NaN; end end end
end6
Output Results to NetCDF File
Output Results to NetCDF File
% Output Results to NetCDF File
% Create NetCDF file and write data
ncid = netcdf.create('E:\北大兼职\臭氧高温复合\extreme_pre\extreme_precipitation_indices_seasonal.nc', 'NETCDF4');
% Define dimensions
dimid_lon = netcdf.defDim(ncid, 'lon', locount);
dimid_lat = netcdf.defDim(ncid, 'lat', lacount);
dimid_season = netcdf.defDim(ncid, 'season', 4);
dimid_year = netcdf.defDim(ncid, 'year', length(Year));
% Define variables
varid_lon = netcdf.defVar(ncid, 'lon', 'double', dimid_lon);
varid_lat = netcdf.defVar(ncid, 'lat', 'double', dimid_lat);
varid_year = netcdf.defVar(ncid, 'year', 'double', dimid_year);
varid_prcptot = netcdf.defVar(ncid, 'PRCPTOT', 'double', [dimid_lon, dimid_lat, dimid_season, dimid_year]);
varid_r10 = netcdf.defVar(ncid, 'R10', 'double', [dimid_lon, dimid_lat, dimid_season, dimid_year]);
varid_r20 = netcdf.defVar(ncid, 'R20', 'double', [dimid_lon, dimid_lat, dimid_season, dimid_year]);
% End definition mode
netcdf.endDef(ncid);
% Write data
netcdf.putVar(ncid, varid_lon, lon);
netcdf.putVar(ncid, varid_lat, lat);
netcdf.putVar(ncid, varid_year, Year);
netcdf.putVar(ncid, varid_prcptot, PRCPTOT_seasonal);
netcdf.putVar(ncid, varid_r10, R10_seasonal);
netcdf.putVar(ncid, varid_r20, R20_seasonal);
% Close NetCDF file
netcdf.close(ncid);
-THE END-
