%% Important notes
% 1) This .m file is a supplement to Deng et al. (2022, Nature Communications). 
%    Manuscript title: A global temperature control of silicate weathering intensity.
% 2) This file provides a template for extracting river basin outlines and basin-scale environmental variables 
%    including topographic metrics and mean annual temperature (MAT) as examples.
% 3) In this file, two variables ('Sample' and 'Loc', see below) can be changed to calculate basin-scale metrics of a target sample. 
% 4) Before running the code below, Matlab code package TopoToolBox 2 should be downloaded 
%    (available at https://github.com/wschwanghart/topotoolbox; Schwanghart & Scherler, 2014) 
% 5) Please read "Readme.docx" before runing the code.

%% Code
clear, clc, close all
% 1. Input variables (can be changed)
Sample = "Example";%Sample name
Loc = [121.2502 28.7533];%[Lon. Lat.]

% 2. Add folder to search path
currentFolder=pwd;% Use the current working directory or specify it by yourself.
% TopoToolbox-related files in folder ...\topotoolbox-master
addpath ([currentFolder '\topotoolbox-master\'])
addpath ([currentFolder '\topotoolbox-master\utilities\'])
addpath ([currentFolder '\topotoolbox-master\IOtools\'])
addpath ([currentFolder '\topotoolbox-master\colormaps\'])
% Grid data of mean annual temperature (rename it as 'MAT.tif') in folder ...\Grid data
addpath ([currentFolder '\Grid data\'])
% 3. Set the path for saving output files
Filepath = ([currentFolder '\Results\']);
warning off
% 4. Load mean annual surface temperature .tif file. Download link: https://www.worldclim.org/data/worldclim21.html
MAT_map = GRIDobj('MAT.tif');
disp('MAT grid data loading finished')
% 5. Download DEM of the study area. Please check internet connection in advance.
% 5.1 Set the extent of DEM for downloading. Note: both numbers can be changed depending on catchment size
Lat_extent = 4; %Unit: degree
Lon_extent = 5; %Unit: degree
% 5.2 Download DEM of the study region
% Please note that the readopentopo tool may need the input of 'apikey' now (not required before). Please type help readopentopo for instructions.
DEM = readopentopo('extent',[round(Loc(1)-Lon_extent,0) round(Loc(1)+Lon_extent,0) round(Loc(2)-Lat_extent,0) round(Loc(2)+Lat_extent,0)],'demtype','SRTMGL3'); %Boundary condition: [West East South North]; DEM resolution: 90 m.
[DEM,utmzone] = reproject2utm(DEM,90);%DEM projection.
DEM.Z(DEM.Z<=0) = nan;
DEM = inpaintnans(DEM);
% 6. Extract the outline of the sampling upstream basin based on the location (Loc).
FD = FLOWobj(DEM,'preprocess','carve');% Flow direction
A = flowacc(FD); % Flow accumulation
DEM = imposemin(FD,DEM,0.0001);%Minima imposition along stream network
S = STREAMobj(FD,'minarea',1e6/A.cellsize^2);% Create stream object; upslope area threshold (1 km2) can be modified. 
[xs,ys] = mfwdtran(DEM.georef.mstruct,Loc(2),Loc(1));% Project sampling location
[xs,ys] = snap2stream(S,xs,ys);% Snap the sampling location to the nearest stream location because the reported location may not be accurate.
% The code above can be modified by adding more constraints (e.g. 'streamorder','>=3'; see 'help STREAMobj/snap2stream'). 
M = drainagebasins(FD,xs,ys);%Obtain the upstream basin from the sampling point
disp('Basin outline extraction finished')
% 7. Calculate topographic metrics
% 7.1 Calculate the centroid of the studied basin
STATS = regionprops(M.Z,'Centroid');
[lat_cen,lon_cen] = ind2coord(M,sub2ind(M.size,round(STATS.Centroid(2)),round(STATS.Centroid(1))));% Unit: m
[lat_cen_deg,lon_cen_deg] = minvtran(DEM.georef.mstruct,lat_cen,lon_cen);% Unit: degree
% 7.2 Calculate flow length (FL)
DEMc = crop(DEM,M,NaN);
FDc = FLOWobj(DEMc,'preprocess','c');
FAc = flowacc(FDc);
Sc = STREAMobj(FDc,'minarea',1e6/FAc.cellsize^2);
Sc = klargestconncomps(Sc,1);
Dc = flowdistance(FDc,'downstream');
Dc = Dc/1000;
FL = nanmax(Dc(:));%Unit: km
% 7.3 Calculate drainage area (Area)
a = FAc;
a.Z = a.Z.*(a.cellsize^2)./10^6;
Area = nanmax(a.Z(:));%Unit: km2
% 7.4 Calculate mean elevation (Meanele)
Ele = DEMc.Z;
Meanele = nanmean(Ele(:));%Unit: m
% 7.5 Calculate mean local slope
Slopeobj = gradient8(DEMc,'degree');
MeanSlope = nanmean(Slopeobj.Z(:));%Unit: deg
% 7.6 Calculate max. relief
Maxrelief = nanmax(Ele(:))-nanmin(Ele(:));% m
% 8. Calculate basin-averaged mean annual surface temperature
MATc = crop(MAT_map,[round(Loc(1)-Lon_extent,0) round(Loc(1)+Lon_extent,0)],[round(Loc(2)-Lat_extent,0) round(Loc(2)+Lat_extent,0)]);
MATc = reproject2utm(MATc,DEM,'method','bilinear');
MAT_clip = crop(MATc,M,NaN);
MAT_FAc = flowacc(FDc,MAT_clip);
MATc_area = MAT_clip;
MATc_area.Z(isnan(MAT_clip.Z)==0) = 1;
MATc_area_FAc = flowacc(FDc,MATc_area);
outlet_ind = (MATc_area_FAc.Z==nanmax(MATc_area_FAc.Z(:)));
MAT_Aweighted = GRIDobj(MAT_clip);
MAT_Aweighted.Z = MAT_FAc.Z./MATc_area_FAc.Z;
MAT_data = MAT_Aweighted.Z(outlet_ind);% Basin-averaged MAT, ^{o}C
% 9. Data export (.xlsx)
Title = {'Sample name' 'Longitude (deg)' 'Latitude (deg)' 'Centroid lon. (deg)' 'Centroid lat. (deg)' 'Flow length (km)' 'Drainage area (km2)' 'Mean elevation (m)' 'Mean local slope (deg)' 'Max. relief (m)' 'MAT (deg)'};
xlswrite([Filepath char(Sample) '.xlsx'],Title,'Result','A1');
MX = [Sample Loc(1) Loc(2) lon_cen_deg lat_cen_deg FL Area Meanele MeanSlope Maxrelief MAT_data];
xlswrite([Filepath char(Sample) '.xlsx'],MX,'Result','A2');
disp('Data saved')
% 10. Check if the extracted basin is correct and export figure (.pdf)
% If the basin outline differs from that reported in literature, please check the accuracy of given location data or add more constraints on 'snap2stream' function 
figure
hold on 
imageschs(DEM,[],'colormap','landcolor')
mask = DEMc;
mask.Z(mask.Z>0) = 9999;
trans = ~isnan(mask);
imagesc(mask,'AlphaData',trans.Z)
plot(Sc,'w','LineWidth',0.2)
[x_original,y_original] = mfwdtran(DEM.georef.mstruct,Loc(2),Loc(1));
scatter(x_original,y_original,'filled','b')
ax = gca;
ax.Box = 'on';
cb = colorbar;
cb.Label.String = 'Elevation (m)';
saveas(gcf,[Filepath char(Sample) '.pdf'])

