chadagreene / cdt Goto Github PK
View Code? Open in Web Editor NEWThe Climate Data Toolbox for MATLAB
Home Page: https://www.chadagreene.com/CDT/CDT_Contents.html
License: Apache License 2.0
The Climate Data Toolbox for MATLAB
Home Page: https://www.chadagreene.com/CDT/CDT_Contents.html
License: Apache License 2.0
The anomaly plot function could be much more versatile by allowing the user to specify a variable baseline. This is useful for visualising differences between the annual trend of a given parameter and its climatological mean. This is an example plot:
Implementing this into the code is pretty straightforward:
Assert that the base is either a scalar (constant baseline) or a vector the same length as x and y (variable baseline):
assert(isscalar(base) | numel(base)==numel(y),'Input error: Base value must be a scalar or a vector of the same length as x and y.')
If base is a scalar, convert it into a column vector the length of y (before "Columnate inputs to ensure...":
% If base is a scalar convert it into a column vector the length of y:
if isscalar(base)
base = repmat(base,numel(y),1);
end
Columnate base just like x and y (only relevant if base is not a scalar):
base = base(:);
Find NaNs both in y and base so filling will work (and remove them from base as well):
% If y or base contains nans, ignore them so filling will work:
ind = (isfinite(y) & isfinite(base));
...
base = base(ind);
The intersections
subfunction now simplifies as follows:
[xc,yc] = intersections(x,y,x,base);
Add zero crossings to base:
% Add zero crossings to the input dataset and sort them into the proper order:
...
base = [base;yc];
...
base = base(ind);
Splitting the data into a top and bottom dataset changes as follows:
% Start thinking about this as two separate datasets which share baseline values where they meet:
ytop = y;
ytop(ytop<base) = base(ytop<base);
ybot = y;
ybot(ybot>base) = base(ybot>base);
Instead of using area
we now need to use fill
, the output remains the same:
% Plot the top half:
htop = fill([x;flipud(x)],[ytop;flipud(base)],topcolor);
% Plot the bottom half:
hbot = fill([x;flipud(x)],[ybot;flipud(base)],bottomcolor);
VOILA! This will allow the user to provide either a constant or a variable baseline and the code will do the rest:
Here is a copy of the adjusted code:
anomaly_JW.m.zip
I hope this is helpful.
Cheers,
Jake
SPEI returns a time step equal to the input variable integrationtime
. Although not wrong, that does not agree with the spirit of SPEI. The time step of the returned SPEI time-series should be the same as the meteorological time-series that goes as input. The integration should then occur as a sliding window.
CDT 中文翻译.pdf
I'm using CDT and I'm exciting about finding such a terrific toolbox for processing meteorology data. However as I'm a Chinese, some of my colleagues or classmates cannot understand the wiki and instruction exactly. So I translate the CDT Contents (mostly from https://www.chadagreene.com/CDT/CDT_Contents.html ) into Chinese.
The work has been started about 5 months ago. Due to my personal reasons, I feel very ashamed of delaying the process of translation. Last week, I translated and typeset all the function documents of the toolbox, which took a long time longer than expected. This is the PDF version after typesetting. I'm not sure if there will be any errors in it. Please review it if convenient.
I hope it may help other people when using this toolbox.
Hi,
I'm new to CDT and new to GitHub. I am trying to write my thesis on climate change with the help of Matlab and CDT…
And I found out that the geomask does not work well when entering multiple countries = polygons (at least according to the instructions).
Therefore, instead of the part in the geomask.m function:
indv = lonv > 180;
lonv(indv) = lonv(indv) - 360;
I suggest using:
% The following part of the code (lonv correction) was made by Tomáš Poláček
% Check if lonv is of type cell
if iscell(lonv)
% Iterate over individual rows
for i = 1:numel(lonv)
% If the given line contains double data
if isnumeric(lonv{i})
indv = lonv{i} > 180;
lonv{i}(indv) = lonv{i}(indv) - 360;
end
end
else
% If lonv is not of cell type, it is a different variable type,
% we perform the operation on the entire lonv (if it is of type double)
if isnumeric(lonv)
indv = lonv > 180;
lonv(indv) = lonv(indv) - 360;
end
end
There is also inaccuracy in CDT help (geomask):
% For this sample 1 degree resolution grid:
[Lat,Lon] = geogrid;
(geogrid is not defined function)
So maybe it will be better to use something like:
[Lon,Lat] = meshgrid(lon, lat);
Thank you in advance for your response and possible correction of my statement.
If coordinates on a 0 to 360
grid cross 180
then the resulting mask will nearly be the inverse of the intended mask. The reason for this is because of the combination of changing to a -180 to 180
coordinate system and the assumption that the smaller longitude is the leftmost and larger rightmost (lines referenced below).
Lines 95 to 99 in 8838903
Line 141 in 8838903
EX (simplified to 1d):
lon = [0 60 120 180 240 300];
lonv = [100 250];
the expected mask is
[0 0 1 1 0 0]
both inputs would be changed to -180 to 180
coordinate system
lon = [0 60 120 180 -120 -60];
lonv = [100 -110];
min(lonv)
is now -110
instead of 100
and lon>min(lonv)
would produce
[1 1 1 1 0 1]
and lon<max(lonv)
would produce
[1 1 0 0 1 1]
the logical and of these 2 produces
[1 1 0 0 0 1]
I propose assuming the first value of lonv is the leftmost input. This would be more in line with how someone would read the value and find it on a globe.
if lonv(1) > lonv(2)
mask = lat>min(latv) & lat<max(latv) & (lon>lonv(1) | lon<lonv(2));
else
mask = lat>min(latv) & lat<max(latv) & (lon>lonv(1) & lon<lonv(2));
end
The current version of daily_insolation.m in the CDT calculates the orbital parameters from the BER90 solution described in Berger & Loutre 1991. As noted in Berger & Loutre 1991 (see conclusion bullet point no. 4) the BER78 solution "is preferable to BER90 for the last glacial-interglacial cycles because of its better accuracy close to present-day times". This is most easily checked by comparing the eccentricity value from daily_insolation.m (0.0172571) for 2000CE, compared to Wikipedia (0.0167086), a difference of 0.5485e-3. Using the BER78 solution eccentricity = 0.0167441, a difference of 0.0355e-3 (compared to Wikipedia) and an extra degree of precision. A similar improvement is found in the longitude of perihelion for 2000CE, Wikipedia = 102.9, daily_insolation.m = 100.51, and BER78 = 101.18.
I therefore propose the attached solution (daily_insolation_new.zip) which requires the mat file Berger.zip. This new function allows the user to pass in an option handle 'BER78' which calculates the eccentricity, obliquity, and longitude of perihelion from the BER78 solution described in Berger & Loutre 1991.
Description of changes
I have updated the comments under 'syntax' and 'description' to describe the new functionality.
The updates to the code are:
line 80 to handle the case of multiple options, i.e. 'constant',value, 'mjmd', and 'BER78'.
narginchk(3,8)
lines 107-111 to handle the new option
if any(strcmpi(varargin,'BER78'))
BER78 = true; % calculate orbital parameters from Berger 1978
else
BER78 = false;
end
if BER78
[ecc,epsilon,omega]=orbital_parameters_BER78(kyear); % function is below in this file
else
[ecc,epsilon,omega]=orbital_parameters(kyear); % function is below in this file
end
% === Calculate orbital parameters from Berger 78 ===
% as described in Berger & Loutre 1991
function [ecc,epsilon,omega] = orbital_parameters_BER78(kyear)
t=-kyear.*1000; % convert kyr to yr and flip sign
% Load the tables from Berger & Loutre 1991
load('Berger.mat','Berger');
Table1 = Berger.Berger78.Table1;
Table4 = Berger.Berger78.Table4;
Table5 = Berger.Berger78.Table5;
% Constants from Table 9 in Berger & Loutre 1991
psibar = 50.439273/60./60. * pi/180 ; % kbar
estar = 23.320556; % eps *
zeta = 3.392506 * pi/180; % alpha
twopi = 2*pi;
sectorad = pi/(180*60.*60);
% Order of table columns 'Term','Amp','Rate','Phase','Period'
M = Table4(:,2); % Amp
g = Table4(:,3).*sectorad; % Rate
b = Table4(:,4).*pi./180; % Phase
F = Table5(:,2).*sectorad; % Amp
fp = Table5(:,3).*sectorad; % Rate
d = Table5(:,4).*pi./180; % Phase
A = Table1(:,2)./60./60; % Amp
f = Table1(:,3).*sectorad; % Rate
phi = Table1(:,4).*pi./180; % Rate
eps = estar + sum(A.*cos(f.*t+phi));
esinpi = sum(M.*sin(g.*t+b));
ecospi = sum(M.*cos(g.*t+b));
psi = psibar.*t + zeta + sum(F.*sin(fp.*t +d));
e = sqrt(esinpi.^2+ecospi.^2);
Pi = atan(esinpi./ecospi)+pi.*(ecospi<0);
eps = eps * pi./180;
varpi = mod(Pi+psi+pi,twopi);
epsilon=eps;
ecc=e;
omega=varpi;
end
In addition to the above, it would not be too difficult to calculate the BER90 solution in the same way as the BER78 method
orbital_parameters_BER90.zip.
I hope you find this useful.
Andrew
User Behzad Navidi asked if implementing SPI would require significant additional changes to the code.
Answer:
Computing the SPI will only require a small
change in the code. The water balance D=P-E will no longer be necessary
and can be replaced with D=P. This would mean that the standardization
of the variable later could be made with the gamma function, but I
don't see a reason why the currently implemented log-logistic won't work.That would require creating a new function spi.m that would be similar
to spei.m only with the small change described above. Definitely worth
doing.
ind = S>=0;
Line 121: Z(ind) = ((S(ind)-1)/StdS).(S(ind));
should be: Z(ind) = ((S(ind)-1)/StdS).(S~=0);
The (S(ind)) is only correct when ind = logical 0. When ind = logical 1, Z(ind) = ((S(ind)-1)/StdS).*(S(ind)) = ((S(ind)-1)/StdS).*S. This is incorrect.
In the function nao
, the mean is removed in lines 64 through 66. However, the mean is already removed by the standardize
call in lines 60 through 62, which makes the second removal redundant.
sw_pden
in mld
calls two non-existent functions sw_ptmp
and sw_dens
.
It would be nice to add an option to ekman.m to use directly provided wind stress rather than using its internal calculation through calling windstress() function. Anyway, this is just a feature request and I think it is fairly easy to implement.
Starting in R2022a, the builddocsearchdb function creates the subfolder helpsearch-v4 to contain the search database files. Previously, builddocsearchdb created a subfolder named helpsearch-v3. As a Result the CDT Documentation is not searchable in Matlab R2022a or later
.
The constant offset (termoption 3 or 4) is not correct.
Hello Chad Greene,
I'm not sure if this page is still active and if you're responding to the queries.
But, I have a few basic related issues:
My inputs are 3*3 matrix: 161 321 66: 66 is the time component !!! (seasonal data - 3 months for 22 years).
Can the team please help me correct this?
Appreciate your work! It's a useful toolbox!
In mld.m profile, IN-SITU temperature is inputted. However, to calculate potential density, the code refers to gsw_rho, in which the input arguments are CONSERVATIVE temperature (or potential temperature) and ABSOLUTE SALINITY (g/kg). The unit of salinity is not illustrated in the header lines. The difference in the type of temperature and an uncommon salinity unit could lead to a serious mistake. Simply change gsw_rho to sw_pden would do the trick. :p
ps: There's a typo on https://www.chadagreene.com/CDT/mld_documentation.html. The examples adopt a T threshold of 0.2 deg C, yet it's changed to 0.02 deg in the code below.
Hi Chad,
I am writing a navigating code for a vessel in seaway, and one of the requirements is to check if the next coordinate is land or sea. If I know a set of coordinates of subsea constructions, oil rigs, etc. How can I add these polygon area into the island function.
Thank you for your time and concern reading my question,
Sicnerely,
The distance2coast.mat file contains the distance of 1/8 degree grid over ocean from the coast. But at the above mentioned location it do not follow. The distance2coast.mat data is verified with load coast data and Etopo2 data. In both the cases it do not follow.
Known issues i noted while creating this issue:
"this function is probably best suited for mesoscale processes that aren't affected by little islands." (Source: https://chadagreene.com/CDT/dist2coast_documentation.html)
The integrated summer insolation should be integrated from 152 day to 243 day (June to August), not select 'Fsw(Few<275)=0'. If I am right, please revise it.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.