Tuesday, September 15, 2026

Matlab: AQHI

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
% Created by LI Xu
% Version 1.0
% September 10, 2026


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy
% http://lixuworld.blogspot.com/

% https://aqihub.info/indices/canada

clear;
clc;
close all;

timebegin=tic;
cur_data=date;
cur_time=fix(clock);
str1=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');

%% ============================================================
%  Canadian AQHI from NO2, O3 and PM2.5 GeoTIFF rasters
%
%  Required units:
%     NO2  : ppb
%     O3   : ppb
%     PM2.5: ug/m^3
%
%  Ideally each input raster is a 3-hour rolling mean.
%% ============================================================



% Source Directory
SouDir='input';

no2files=dir(fullfile(SouDir, 'NO2*'));
o3files=dir(fullfile(SouDir, 'O3*'));
pm25files=dir(fullfile(SouDir, 'PM*'));


for ii=1:numel(no2files)

    no2path=fullfile(SouDir, no2files(ii).name);
    o3path=fullfile(SouDir, o3files(ii).name);
    pm25path=fullfile(SouDir, pm25files(ii).name);

    [~, fname, ~]=fileparts(pm25path);
    fname=strsplit(fname, '_');
    fname=fname{end};

    otpath=fullfile(SouDir, ['AQHI_', fname, '.tif']);

    GenAQHICompute(no2path, o3path, pm25path, otpath);


    disp(['[', num2str(ii), '/', num2str(numel(no2files)), ']~', otpath]);





end


fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');
cur_data=date;
cur_time=fix(clock);
str2=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time End: ');
disp(str2);
timespan=toc(timebegin);
fprintf('Time Span: %.4f s\n', timespan);


disp('******************************************************************');

function GenAQHICompute(no2path, o3path, pm25path, otpath)
%% 1. Read pollutant rasters

    [NO2, R_NO2] = readgeoraster(no2path, ...
        'OutputType', 'double');
    
    [O3, R_O3] = readgeoraster(o3path, ...
        'OutputType', 'double');
    
    [PM25, R_PM25] = readgeoraster(pm25path, ...
        'OutputType', 'double');

    %% 2. Make sure the rasters have the same dimensions

    if ~isequal(size(NO2), size(O3), size(PM25))
        error('NO2, O3 and PM2.5 rasters must have the same dimensions.');
    end

    %% 3. Handle physically impossible negative concentrations
    % Optional but usually sensible for interpolated/modelled concentration maps.

    NO2(NO2 < 0)   = 0;
    O3(O3 < 0)     = 0;
    PM25(PM25 < 0) = 0;

    %% 4. Identify pixels where all three pollutants are available

    valid = isfinite(NO2) & isfinite(O3) & isfinite(PM25);

    AQHI_raw = nan(size(NO2));

    AQHI_raw(valid) = (10 / 10.4) .* 100 .* ...
        ( ...
        (exp(0.000537 .* O3(valid))   - 1) + ...
        (exp(0.000871 .* NO2(valid))  - 1) + ...
        (exp(0.000487 .* PM25(valid)) - 1) ...
        );


    %% 6. Calculate the reported integer AQHI

    AQHI = nan(size(AQHI_raw));

    AQHI(valid) = round(AQHI_raw(valid));


    info = geotiffinfo(no2path);
    geotiffwrite( ...
    otpath, ...
    single(AQHI), ...
    R_NO2, ...
    'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);


end

References

  1. https://www.canada.ca/en/environment-climate-change/services/air-quality-health-index/about.html

Tuesday, August 18, 2026

Video: Air Parcel Trajectory

Two-week HYSPLIT simulations of air parcel trajectories originating in Manitoba and transported across the globe, using ERA5 reanalysis data as meteorological input.

Video: Dust Storm

NOAA GOES-19 imagery showing a dust storm affecting the entire Winnipeg area on the afternoon of May 14, 2026.

Reference:

https://classic107.com/articles/rare-dust-storm-warning-issued-across-southern-manitoba

Thursday, August 6, 2026

Matlab: An example Map script

% Created by LI Xu
% Version 1.0
% August 6, 2026


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy
% http://lixuworld.blogspot.com/


clear;
clc;
close all;

timebegin=tic;
cur_data=date;
cur_time=fix(clock);
str1=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');



% Define the Google Street Map tile URL
% name = 'google_streets';
name='streets';
% url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/${z}/${y}/${x}';
% url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/${z}/${y}/${x}';
% url = 'https://mt1.google.com/vt/lyrs=m&x=${x}&y=${y}&z=${z}&style=p.v:off|s.t:3|s.e:l|p.v:on';
% 
% url = 'https://mt1.google.com/vt/lyrs=m&x=${x}&y=${y}&z=${z}&style=p.v:off|s.t:3|s.e:l|p.v:on';
% url = 'https://mt1.google.com/vt/lyrs=m&x=${x}&y=${y}&z=${z}';
% url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Physical_Map/MapServer/tile/${z}/${y}/${x}';
% attribution = 'Map data ©2026 Google';
% url = 'https://cartodb-basemaps-a.global.ssl.fastly.net/light_all/${z}/${x}/${y}.png';
% url = 'https://cartodb-basemaps-a.global.ssl.fastly.net/dark_all/${z}/${x}/${y}.png';
% name = 'carto_voyager';
% attribution = '© OpenStreetMap contributors, © CARTO';
% url = 'https://cartodb-basemaps-a.global.ssl.fastly.net/rastertiles/voyager/${z}/${x}/${y}.png';
% 
% 
% % Add it to your MATLAB basemap list
% addCustomBasemap(name, url, 'Attribution', attribution);

% url = 'https://server.arcgisonline.com/ArcGIS/rest/services/NatGeo_World_Map/MapServer/tile/${z}/${y}/${x}';
% url = 'https://server.arcgisonline.com/ArcGIS/rest/services/Specialty/NatGeo_World_Map/MapServer/tile/${z}/${y}/${x}';
% url = 'https://cartodb-basemaps-a.global.ssl.fastly.net/rastertiles/voyager/${z}/${x}/${y}.png';
attribution = '© OpenStreetMap contributors, © CARTO';

url = 'https://mt1.google.com/vt/lyrs=s&x=${x}&y=${y}&z=${z}';
addCustomBasemap(name, url, 'Attribution', attribution);
% addCustomBasemap(name, url, 'Attribution', 'Sources: Esri, National Geographic, Garmin');
% url = 'https://stamen-tiles.a.ssl.fastly.net/toner/${z}/${x}/${y}.png';
% addCustomBasemap('stamen_toner', url, 'Attribution', 'Map tiles by Stamen Design, CC BY 3.0');

% Create geographic axes with OpenStreetMap basemap
figure;
gx = geoaxes('Basemap', name);
set(gcf, 'OuterPosition', get(0, 'ScreenSize'));
f = gcf;

% Define custom latitude/longitude range (Winnipeg-centered)
% latRange = [49.76, 50.01];
% lonRange = [-97.45, -96.9];
latRange = [49.79+0.005+0.009, 49.82-0.0068];
lonRange = [-97.15+0.008, -97.13+0.004];
% Calculate aspect ratio for exact limits without expansion
mean_lat = mean(latRange);
cos_factor = cosd(mean_lat);
ar = diff(latRange) / (diff(lonRange) * cos_factor);  % Required height/width for axes

% Set up resize function to maintain aspect and exact limits
gx.Units = 'pixels';
f.Units = 'pixels';
f.SizeChangedFcn = @(~, ~) myresizefcn(gx, f, ar, latRange, lonRange);
myresizefcn(gx, f, ar, latRange, lonRange);
% Rotate the whole map figure visually
view(gx, 2);
% rotate3d(f, 'on');
% % rotate3d(gx, 'on');
% gx.View = [30 90];


gx.Grid = 'off';
gx.TickDir = 'none';
gx.LatitudeLabel.String = '';
gx.LongitudeLabel.String = '';
gx.FontSize = 0.01;


% dat=cell2mat(dat);
% lats=dat(:, 1);
% lons=dat(:, 2);
% nums=dat(:, 3);
% hold on;
% for i = 1:length(lats)
%     % Draw Circle
%     plot(gx, lats(i), lons(i), 'bo', 'MarkerSize', 18, 'MarkerFaceColor', 'y');
%     % Draw Number
%     text(gx, lats(i), lons(i), num2str(nums(i)), ...
%         'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'Color', 'blue', 'FontWeight','bold', 'FontSize', 15);
% end

% 49.81155561273022, -97.13142673818446
plot(gx, 49.81155561273022, -97.13142673818446, 'pentagram', 'MarkerSize', 20, 'MarkerFaceColor', 'blue');


hold off;

gx.Grid = 'off';
gx.TickDir = 'none';
gx.LatitudeLabel.String = '';
gx.LongitudeLabel.String = '';
gx.FontSize = 0.01;

% dim = [0.1 0.83 0.8 0.1]; % Position at the top center
% str = 'Campus Water Fountain Map';
% 
% hAnn = annotation('textbox', dim, 'String', str, ...
%     'EdgeColor', 'none', ...         % No border
%     'HorizontalAlignment', 'center', ...
%     'VerticalAlignment', 'middle', ...
%     'FontSize', 25, ...
%     'FontWeight', 'bold', ...
%     'Color', 'white', ...
%     'Interpreter', 'none');

pause(3);
otname='kunakmap_vector';
export_fig(otname, '-png', '-painters', '-r300');
close all;



fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');
cur_data=date;
cur_time=fix(clock);
str2=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time End: ');
disp(str2);
timespan=toc(timebegin);
fprintf('Time Span: %.4f s\n', timespan);



disp('************************************************************');

function myresizefcn(gx, f, ar, latlim, lonlim)
    if (f.Position(4) / f.Position(3)) < ar
        % Figure is wide: constrain by height
        gx.Position(4) = f.Position(4) * 0.8;
        gx.Position(3) = gx.Position(4) / ar;
        % Center
        gx.Position(1) = (f.Position(3) - gx.Position(3)) / 2 + f.Position(1);
        gx.Position(2) = (f.Position(4) - gx.Position(4)) / 2 + f.Position(2);
    else
        % Figure is tall: constrain by width
        gx.Position(3) = f.Position(3) * 0.8;
        gx.Position(4) = gx.Position(3) * ar;
        % Center
        gx.Position(1) = (f.Position(3) - gx.Position(3)) / 2 + f.Position(1);
        gx.Position(2) = (f.Position(4) - gx.Position(4)) / 2 + f.Position(2);
    end
    % Re-apply limits to maintain exact range
    geolimits(gx, latlim, lonlim);
end

Tuesday, May 12, 2026

Matlab: Table on the PPT created via Matlab

% Created by LI Xu
% Version 1.0
% August 28, 2025


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy
% http://lixuworld.blogspot.com/

clear;
clc;
close all;

timebegin=tic;
cur_data=date;
cur_time=fix(clock);
str1=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');

import mlreportgen.ppt.*

% 1. Create a presentation
ppt = Presentation('AQ.pptx');

% 2. Create a slide (e.g., 'Blank' or 'Title and Content')
slide = add(ppt, 'Title and Content');
replace(slide, 'Title', 'Measured Data Table');

% 3. Prepare your data (Cell array: Rows x Cols)
data = {'Variable', 'Campus', 'Winnipeg'; 'Temperature (℃) ', '40', '15'; ...
    'Relative Humidity (%)', '65', '70'; ...
    'Wind speed (m/s)', '65', '70'; ...
    'PM2.5 (μg/m³)', '65', '70'; ...
    'PM10 (μg/m³)', '65', '70'; ...
    'NO₂ (ppb)', '65', '70'; ...
    'O₃ (ppb)', '65', '70'; ...
    'CO₂ (ppm)', '65', '70'};


% 4. Create and add the table
% 'Table' automatically creates rows/columns based on the data
tableObj = Table(data);
firstRow = tableObj.Children(1);

firstRow.Style = {BackgroundColor('#296080'), FontColor('white'), Bold()};
tableObj.Style = {HAlign('center'), VAlign('middle'), Border('solid'), RowHeight('0.3in')};
% tableObj.ColWidths = {'1.5in', '1.0in', '1.0in'};
% Column 3: 1.0 inches (Winnipeg)
spec1 = ColSpec('2.5 in');
spec2 = ColSpec('1.5 in');
spec3 = ColSpec('1.5 in');

% 2. Assign the array of specs to the table
tableObj.ColSpecs = [spec1, spec2, spec3];


% 1. Target the specific cell (Row 2, Column 2) using entry(row, col)
targetCell = entry(tableObj, 2, 2);

% 2. Apply the background color
targetCell.Style = {BackgroundColor('red')};

% 1. Target the specific cell (e.g., Row 2, Column 2)
targetCell = entry(tableObj, 2, 2);

% 2. Apply the font color (e.g., red)
% You can use names like 'red', 'white', or hex codes like '#FF0000'
targetCell.Style = [targetCell.Style, {FontColor('white')}];


% Add the table to the content placeholder
replace(slide, 'Content', tableObj);

% 5. Close and View
close(ppt);
rptview(ppt);



fprintf('Time Begin: ');
fprintf(str1);
fprintf('\n');
cur_data=date;
cur_time=fix(clock);
str2=sprintf('%s %.2d:%.2d:%.2d', cur_data, cur_time(4), cur_time(5), cur_time(6)); 
fprintf('Time End: ');
disp(str2);
timespan=toc(timebegin);
fprintf('Time Span: %.4f s\n', timespan);

disp('*************************************************************************');