-
Notifications
You must be signed in to change notification settings - Fork 31
/
equistat.m
95 lines (78 loc) · 3.27 KB
/
equistat.m
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
function [lonlon,latlat,intlon,intlat,orlon,orlat]=equistat(c11,cmn,lonnum,latnum)
% [lonlon,latlat,intlon,intlat,orlon,orlat]=EQUISTAT(c11,cmn,lonnum,latnum)
%
% Expresses a geographical grid of regularly spaced latitudes and
% longitudes in terms of the true distances along the surface of the
% sphere. Suitable to interpolate a geographical grid to a locally
% regular Cartesian flat grid.
%
% INPUT:
%
% C11 [lon lat] of the upper left corner of the map [degrees]
% CMN [lon lat] of the lower right corner of the map [degrees]
% lonnum Number of samples across (E-W)
% latnum Number of samples down (N-S)
%
% OUTPUT:
%
% lonlon The Cartesion x-coordinates of the geographical [lon lat] grid
% latlat The Cartesion y-coordinates of the geographical [lon lat] grid
% intlon Regular Cartesian longitudes for interpolation using GRIDDATA
% intlat Regular Cartesian latitudes for interpolation using GRIDDATA
% All of the above are with respect to a top left corner [0,0]
% orlon The longitudes of the original geographical grid
% orlat The latitudes of the original geographical grid
% The above two with respect to the actual top left corner
%
% EXAMPLE I:
%
% XIM=[115 155]; YIM=[-5 -50];
% [lonlon,latlat,intlon,intlat,orlon,orlat]=...
% equistat([XIM(1) YIM(1)],[XIM(2) YIM(2)],20,20);
% fridplot(lonlon,latlat); hold on
% co=fridplot(intlon,intlat); set(co,'Color','b') ; hold off
%
% EXAMPLE II:
%
% XIM=[100 160]; YIM=[1 -50];
%% Get the data on a regular geographical grid (longitude and latitude)
% z=flipud(etopo(fullfile(getenv('IFILES'),'TOPOGRAPHY','EARTH'),5,sort(YIM),XIM));
% subplot(121)
% imagef([XIM(1) YIM(1)],[XIM(2) YIM(2)],z); axis image;
%% Plot the continental outlines, in geographical coordinates, on top
% [a,b,XY]=plotcont([XIM(1) YIM(1)],[XIM(2) YIM(2)]); axis tight
% cb=cax2dem([-8000 1500]); delete(cb)
%% Find interpolation points on a regular Cartesian grid (E-W and N-S)
%% contained in the original grid
% [lonlon,latlat,intlon,intlat]=...
% equistat([XIM(1) YIM(1)],[XIM(2) YIM(2)],size(z,2),size(z,1));
%% Interpolate the geographical grid to this regular Cartesian grid
% zi=griddata(lonlon,latlat,z,intlon,intlat);
% subplot(122)
% Plot the interpolated image on the regular Cartesian grid
% imagef([0 0],[range(intlon(:)) -range(intlat(:))],zi); axis image
%% Project the continental outlines according to the same scheme
% [flon,tlat]=project(XY(:,1),XY(:,2),[XIM(1) YIM(1)],[XIM(2) YIM(2)]);
% cb=cax2dem([-8000 1500]); delete(cb)
% The next line needs to be fixed !!!!
%% hold on ; d=plot(flon...,tlat...,'k'); hold off
%
% See also UNPROJECT, PROJECT, POLARGRID, ELL2CAR, MASTERSET
%
% Last modified by fjsimons-at-alum.mit.edu, 09/24/2008
XIM=[c11(1) cmn(1)];
YIM=[c11(2) cmn(2)];
longrid=linspace(XIM(1),XIM(2),lonnum);
latgrid=linspace(YIM(1),YIM(2),latnum);
[orlon,orlat]=meshgrid(longrid,latgrid);
lonlon=(orlon-c11(1)).*cos(orlat*pi/180);
latlat=(orlat-c11(2));
% Origin in middle of data set
shift=repmat((lonlon(1,lonnum)-lonlon(:,lonnum))/2,1,lonnum);
lonlon=lonlon+shift;
latlat=latlat;
intval=max(lonlon(:,lonnum)-lonlon(:,lonnum-1));
minval=max(lonlon(:,1));
maxval=min(lonlon(:,lonnum));
intlon=repmat(minval:intval:maxval,latnum,1);
intlat=repmat(latlat(:,1),1,size(intlon,2));