-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
30 changed files
with
3,667 additions
and
398 deletions.
There are no files selected for viewing
Large diffs are not rendered by default.
Oops, something went wrong.
Binary file not shown.
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,67 @@ | ||
| % This script defines some common covariance functions. | ||
| % | ||
| % kMat1 - Matern 1/2 | ||
| % kMat3 - Matern 3/2 | ||
| % kMat5 - Matern 5/2 | ||
| % kSE - squared exponential | ||
| % kDP - dot product | ||
| % kCONST - constant | ||
| % kSIN - square exponential of sine function | ||
| % kMODSIN - product of kSIN and kSE | ||
| % kFREQ - cosine function | ||
| % kRQ - rational quadratic | ||
| % kMatG - Matern generic | ||
| % kDELTAG - delta function on a sphere | ||
| % kGEOG - Matern 5/2 on a sphere | ||
| % kGEOGG - Matern generic on a sphere | ||
| % | ||
| % It also defines first and second derivatives kMat3d, kMat3dd and kDPd, kDPdd | ||
| % | ||
| % Last updated by Robert Kopp, robert-dot-kopp-at-rutgers-dot-edu, Sat Nov 28 16:09:12 EST 2015 | ||
|
|
||
| % we will have a vector of observations | ||
| % and a vector of their uncertainty | ||
|
|
||
| % define covariance functions we will use | ||
|
|
||
| defval('refyear',1970); | ||
|
|
||
| dYears = @(years1,years2) abs(bsxfun(@minus,years1',years2)); | ||
| dYears0 = @(years1,years2) (bsxfun(@minus,years1',years2)); | ||
|
|
||
| angd = @(Lat0,Long0,lat,long) (180/pi)*(atan2( sqrt( (cosd(lat) .* sind(long-Long0)).^2 + (cosd(Lat0) .* sind(lat) - sind(Lat0) .* cosd(lat) .* cosd(long-Long0)).^2),(sind(Lat0) .* sind(lat) + cosd(Lat0) .* cosd(lat) .* cosd(long-Long0)))); | ||
|
|
||
| dDist = @(x1,x2) angd(repmat(x1(:,1),1,size(x2,1)),repmat(x1(:,2),1,size(x2,1)),repmat(x2(:,1)',size(x1,1),1),repmat(x2(:,2)',size(x1,1),1))' + 1e6*(bsxfun(@plus,x1(:,1)',x2(:,1))>1000); | ||
|
|
||
| kMat1 = @(dx,thetas) thetas(1).^2 .* (1).*exp(-dx/thetas(2)); | ||
|
|
||
| kMat3 = @(dx,thetas) thetas(1).^2 .* (1 + sqrt(3)*dx/thetas(2)).*exp(-sqrt(3)*dx/thetas(2)); | ||
|
|
||
| kMat5 = @(dx,thetas) thetas(1).^2 .* (1 + (sqrt(5)*dx/thetas(2)).*(1 + sqrt(5)*dx/thetas(2)/3)).*exp(-sqrt(5)*dx/thetas(2)); | ||
|
|
||
| kSE = @(dx,thetas) thetas(1).^2 * exp(-(dx.^2)/(2*thetas(2).^2)); | ||
|
|
||
| kDELTA = @(dx,thetas) thetas(1).^2 .* (dx==0); | ||
|
|
||
| kDP = @(years1,years2,thetas) thetas(1).^2 * bsxfun(@times,(years1-refyear)',(years2-refyear)); | ||
| kCONST = @(thetas) thetas(1).^2; | ||
| kSIN = @(dt,thetas) thetas(1).^2 * exp(-2*sin(pi*dt/(thetas(2))).^2/thetas(3).^2); | ||
| kMODSIN = @(dt,thetas) kSIN(dt,thetas(1:3)) .* kSE(dt,[1 thetas(4)*thetas(2)]); | ||
| kFREQ = @(dt,thetas) thetas(1).^2 * cos(2*pi*dt/thetas(2)); | ||
| kRQ = @(dx,thetas) thetas(1).^2 * (1 + dx.^2/(2*thetas(2)*thetas(3))).^-thetas(3); | ||
| kMatG = @(dx,thetas) thetas(1).^2 .* 2.^(1-thetas(3))./gamma(thetas(3)) .* (sqrt(2*thetas(3))*(dx+eps)/thetas(2)).^thetas(3) .* besselk(thetas(3),sqrt(2*thetas(3))*(dx+eps)/thetas(2)); | ||
|
|
||
| kDELTAG = @(ad,thetas)thetas(1).^2.*(abs(ad)<1e-4).*(ad<360); | ||
|
|
||
| kGEOG = @(ad,thetas) kMat5(ad,thetas) .* (ad<360); | ||
| kGEOGG = @(ad,thetas)kMatG(ad,thetas).*(ad<360); | ||
|
|
||
|
|
||
| kMat3d = @(years1,years2,dx,thetas)thetas(1).^2.*(-3/(thetas(2).^2)).*dx.*exp(-sqrt(3)*dx/thetas(2)).*(-1+2*bsxfun(@ge,years1',years2)); | ||
| kDPd = @(years1,years2,thetas) thetas(1).^2 * repmat((years1-refyear)',length(years2),1); | ||
| kMat3dd = @(dx,thetas)thetas(1).^2.*(3/(thetas(2).^2)).*(1-sqrt(3)*dx/thetas(2)).*exp(-sqrt(3)*dx/thetas(2)); | ||
| kDPdd = @(years1,years2,thetas) thetas(1).^2 * ones(length(years2),length(years1)); | ||
|
|
||
|
|
||
|
|
||
| FiniteMask = @(x1,x2) (bsxfun(@plus,sum(abs(x1),2)',sum(abs(x2),2)))<1e12; |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,105 @@ | ||
| function [TGdata2,TGdata,thetL,TGmodellocal] = GPSmoothNearbyTideGauges(targcoords,addlsites,thetL0,winlength,thinlength,minlength,optimizemode,psmsldir,gslfile,maxdist,noiseMask) | ||
|
|
||
| % [TGdata2,TGdata,thetL,TGmodellocal] = GPSmoothNearbyTideGauges(targcoords, | ||
| % [addlsites],[thetL0],[winlength],[thinlength],[minlength],[optimizemode], | ||
| % psmsldir,gslfile,[maxdist],[noiseMask]) | ||
| % | ||
| % Find tide gauge sites to include, based on length and proximity criteria, then | ||
| % fit GP model to them in order to interpolate and take running averge. | ||
| % | ||
| % Specifically, from the site of all tide gauge sites with length of at least | ||
| % minlength(3), identifies sites that fit one of the following criteria: | ||
| % 1. are longer than minlength(1) | ||
| % 2. are within maxdist degrees of one of the sites listed in targcoords | ||
| % and are longer than minlength(2) | ||
| % 3. is the nearest tide gauge to one of the sites listed in targcoords | ||
| % 4. has a PSMSL ID listed in addlsites | ||
| % | ||
| % Data are then smoothed and thinned by | ||
| % | ||
| % [TGdata2,thetL,TGmodellocal] = GPSmoothTideGauges(TGdata,winlength, ... | ||
| % optimizemode,thinlength,thetL0); | ||
| % | ||
| % EXAMPLE: | ||
| % | ||
| % [TG,TG0,thetL,TGmodellocal] = GPSmoothNearbyTideGauges(PX.sitecoords, ... | ||
| % [],[],[],[],[150 75 20],1.0,psmsldir,'none'); | ||
| % | ||
| % Last updated by Robert Kopp, robert-dot-kopp-at-rutgers-dot-edu, Thu Nov 06 21:04:21 EST 2014 | ||
| % | ||
|
|
||
| defval('minlength',[150 75 20]); % minimum length for global curves, most purposes and minimum | ||
| defval('thetL0',[]); | ||
| defval('winlength',11); | ||
| defval('thinlength',winlength-1); | ||
| defval('optimizemode',1.1); % set to 1.0 for local optimization only | ||
| defval('maxdist',5); | ||
| defval('thinyrstart',1700); | ||
| defval('noiseMask',[]); | ||
|
|
||
| defval('targcoords',[ | ||
| 34.97 -76.38; %Tump Point, NC | ||
| 35.89 -75.64; %Sand Point, NC | ||
| 39.085 -74.811; % Cape May Courthouse, NJ | ||
| 39.50 -74.415; % Leeds Point, NJ | ||
| 30.6 -81.7; % Nassau, FL | ||
| 44.72 -63.23; % Chezzetcook, Nova Scotia | ||
| ]); | ||
|
|
||
| defval('addlsites',[]); | ||
|
|
||
| defval('psmsldir',fullfile('.','rlr_annual')); | ||
| defval('gslfile',fullfile('.','CSIRO_Recons_gmsl_yr_2011.csv')); | ||
|
|
||
| % read PSMSL data for all sites longer than minlength(3) | ||
| [TGcoords,TGrsl,TGrslunc,TGid,TGsiteid,sitenames,TGsitecoords,sitelen]=ReadPSMSLData(0,1000,minlength(3),psmsldir,gslfile); | ||
|
|
||
| % find long tide gauge sites anywhere | ||
| sub1=find((sitelen>minlength(1))); | ||
| sub1=union(sub1,find(TGsiteid==0)); | ||
|
|
||
| % compute angular distances between sites, and add any tide gauges within | ||
| % maxdist degrees with length greater than minlength(2) | ||
|
|
||
| angd= @(Lat0,Long0,lat,long) (180/pi)*(atan2(sqrt((cosd(lat).*sind(long-Long0)).^2+(cosd(Lat0).*sind(lat)-sind(Lat0).*cosd(lat).*cosd(long-Long0)).^2),(sind(Lat0).*sind(lat)+cosd(Lat0).*cosd(lat).*cosd(long-Long0)))); | ||
|
|
||
| dDist=@(x1,x2)angd(repmat(x1(:,1),1,size(x2,1)),repmat(mod(x1(:,2),360),1,size(x2,1)),repmat(x2(:,1)',size(x1,1),1),repmat(mod(x2(:,2),360)',size(x1,1),1))'+1e6*(bsxfun(@plus,x1(:,1)',x2(:,1))>1000); | ||
|
|
||
| sub2=[]; | ||
| for ii=1:size(targcoords,1) | ||
| TGdist=dDist(TGsitecoords,targcoords(ii,:)); | ||
| [m,mi]=min(TGdist); | ||
| mi=union(mi,find((TGdist<maxdist).*(sitelen'>minlength(2)))); | ||
| sub2=union(sub2,mi); | ||
| end | ||
|
|
||
|
|
||
| % add addlsites | ||
| addlsites=addlsites(:)'; | ||
| subaddlsites=find(ismember(TGsiteid,addlsites)); | ||
|
|
||
| sitesub=union(sub1,sub2); | ||
| sitesub=union(sitesub,subaddlsites); | ||
| sub=find(ismember(TGid,TGsiteid(sitesub))); | ||
|
|
||
| clear TGdata; | ||
| TGdata.datid=TGid(sub); | ||
| TGdata.time1=TGcoords(sub,3); | ||
| TGdata.time2=TGcoords(sub,3); | ||
| TGdata.meantime=TGcoords(sub,3); | ||
| TGdata.limiting=zeros(size(TGid(sub))); | ||
| TGdata.Y=TGrsl(sub); | ||
| TGdata.dY =TGrslunc(sub); | ||
| TGdata.compactcorr=zeros(size(TGid(sub))); | ||
| TGdata.istg = ones(size(TGid(sub))); | ||
| TGdata.lat=TGcoords(sub,1); | ||
| TGdata.long=TGcoords(sub,2);; | ||
| TGdata.Ycv=sparse(diag(TGrslunc(sub).^2)); | ||
| TGdata.siteid=TGsiteid(sitesub); | ||
| TGdata.sitenames=sitenames(sitesub); | ||
| TGdata.sitecoords=TGsitecoords(sitesub,:); | ||
| TGdata.sitelen=sitelen(sitesub); | ||
|
|
||
| %%%%%%%%%% | ||
|
|
||
| [TGdata2,thetL,TGmodellocal] = GPSmoothTideGauges(TGdata,winlength,optimizemode,thinlength,thetL0,thinyrstart,noiseMask); |
Oops, something went wrong.