source: issm/trunk/src/m/classes/boundary.m@ 25836

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

merged trunk-jpl and trunk for revision 25834
File size: 7.1 KB
RevLine 
[24001]1%BOUNDARY class definition
2%
3% Usage:
4% boundary=boundary();
5
6classdef boundary
7 properties (SetAccess=public)
8 shppath = '';
9 shpfilename = '';
10 orientation = 'normal'; %other value is 'reverse'
[25836]11 proj = epsg2proj(4326); %Proj4.0 projection string. , default is WGS 84 (corresponds to epsg 4326)
12
[24001]13 end
14 methods (Static)
15 function self = loadobj(self) % {{{
16 % This function is directly called by matlab when a model object is
17 % loaded. Update old properties here
18 end% }}}
19 end
20 methods
21 function self = boundary(varargin) % {{{
[25836]22 self=setdefaultparameters(self);
23
24 if nargin > 0,
25 options=pairoptions(varargin{:});
26
27 %recover field values:
28 self.shppath=getfieldvalue(options,'shppath','');
29 self.shpfilename=getfieldvalue(options,'shpfilename','');
30 self.orientation=getfieldvalue(options,'orientation','normal');
31
32 %figure out projection string:
33 if exist(options,'epsg'),
34 if exist(options,'proj'),
35 error(['Error in constructor for boundary ' self.shppath '. Cannot supply epsg and proj at the same time!']);
36 end
37 epsg=getfieldvalue(options,'epsg');
38 proj=epsg2proj(epsg); % convert to PROJ.4 format
39 elseif exist(options,'proj'),
40 if exist(options,'epsg'),
41 error(['Error in constructor for boundary ' self.shppath '. Cannot supply proj and epsg at the same time!']);
42 end
43 proj=getfieldvalue(options,'proj');
44 else
45 proj= '+proj=longlat +datum=WGS84 +no_defs';
46 end
47 self.proj=proj;
[24001]48 end
49 end % }}}
50 function self = setdefaultparameters(self) % {{{
51 self.shppath='';
52 self.shpfilename='';
53 self.orientation='normal';
[25836]54 self.proj='+proj=longlat +datum=WGS84 +no_defs'; %EPSG 4326
[24001]55 end % }}}
56 function disp(self) % {{{
57 disp(sprintf(' boundary parameters:'));
58
59 fielddisplay(self,'shppath','path to filename for this boundary');
60 fielddisplay(self,'shpfilename','shape file name');
61 fielddisplay(self,'orientation','orientation (default is ''normal'', can be ''reverse'')');
[25836]62 fielddisplay(self,'proj','shape file projection string (follows PROJ.4 standard)');
[24001]63
64 end % }}}
65 function output=name(self) % {{{
66 output=self.shpfilename;
67 end % }}}
68 function output=edges(self) % {{{
69
[25836]70 %read domain:
71 [path,name,ext]=fileparts(self.shpfilename);
72 if ~strcmpi(ext,'shp'),
73 ext='shp';
74 end
75 output=shpread([self.shppath '/' name '.' ext]);
76
77 %do we reverse?
78 if strcmpi(self.orientation,'reverse'),
79 output.x=flipud(output.x);
80 output.y=flipud(output.y);
81 end
[24001]82 end % }}}
83 function plot(self,varargin) % {{{
84 %recover options
85
86 options=pairoptions(varargin{:});
87
88 %parse input:
89 figurenumber=getfieldvalue(options,'figure',1);
90 color=getfieldvalue(options,'color','r');
91 linewidth=getfieldvalue(options,'linewidth',1);
92 markersize=getfieldvalue(options,'markersize',1);
93 unitmultiplier=getfieldvalue(options,'unit',1);
94 radius=getfieldvalue(options,'radius',6371012);
95 aboveground=getfieldvalue(options,'aboveground',1000);
96 offset=getfieldvalue(options,'offset',.1);
97 fontsize=getfieldvalue(options,'fontsize',10);
[24686]98 label=getfieldvalue(options,'label','on');
[24001]99
[25836]100 if exist(options,'epsg'),
101 if exist(options,'proj'),
102 error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
103 end
104 proj=epsg2proj(getfieldvalue(options,'epsg'));
105 elseif exist(options,'proj'),
106 proj=getfieldvalue(options,'proj');
107 else
108 proj=epsg2proj(4326);
109 end
110
[24001]111 %read domain:
112 [path,name,ext]=fileparts(self.shpfilename);
113 if ~strcmpi(ext,'shp'),
114 ext='shp';
115 end
116 domain=shpread([self.shppath '/' name '.' ext]);
117
[25836]118 %convert boundary to another reference frame: {{{
[24001]119 for i=1:length(domain),
120 try,
[25836]121 [x,y] = gdaltransform(domain(i).x,domain(i).y,self.proj,proj);
[24001]122 catch me
123 disp(me.message());
124 self.disp();
125 end
126 domain(i).x=x; domain(i).y=y;
127 end
128
129 for i=1:length(domain),
130 hold on;
[24686]131 x=domain(i).x*unitmultiplier;
132 y=domain(i).y*unitmultiplier;
[24001]133 if length(x)==1,
134 p=plot(x,y,'k*');
[24686]135 set(p,'MarkerSize',markersize,'Color',color);
136 if strcmpi(label,'on'),
137 t=text(x,y,self.shpfilename,'FontSize',fontsize);
138 end
[24001]139 else
140 p=plot(x,y,'k-');
[24686]141 set(p,'MarkerSize',markersize,'Color',color);
142 if strcmpi(label,'on'),
143 text(sum(x)/length(x),sum(y)/length(y),self.shpfilename,'FontSize',fontsize);
144 end
[24001]145 end
146 set(p,'Color',color);
147 set(p,'LineWidth',linewidth);
148 end
149 %}}}
150 end % }}}
[25836]151 function checkconsistency(self,varargin) % {{{
152 %recover options
153 options=pairoptions(varargin{:});
154 tolerance=getfieldvalue(options,'tolerance',1e-5);
155
156 %recover edges:
157 edges=self.edges();
158
159 if strcmpi(edges.Geometry,'Point'),
160 return;
161 else
162 %check we don't have two identical vertices:
163 x=edges.x; y=edges.y;
164 distance=sqrt( (x(1:end-1)-x(2:end)).^2+(y(1:end-1)-y(2:end)).^2);
165 dmax=max(distance);
166 isidentical=find(distance<dmax*tolerance);
167 if ~isempty(isidentical),
168 error(['boundary ' shpfilename ' has two vertices extremely close to one another']);
169 end
170 end
171 end % }}}
[24001]172 function plot3d(self,varargin) % {{{
173 %recover options
174
175 options=pairoptions(varargin{:});
176
177 %parse input:
178 figurenumber=getfieldvalue(options,'figure',1);
179 color=getfieldvalue(options,'color','r');
180 linewidth=getfieldvalue(options,'linewidth',1);
181 markersize=getfieldvalue(options,'markersize',1);
182 unitmultiplier=getfieldvalue(options,'unit',1);
183 radius=getfieldvalue(options,'radius',6371012);
184 aboveground=getfieldvalue(options,'aboveground',1000);
185 offset=getfieldvalue(options,'offset',.1);
186 fontsize=getfieldvalue(options,'fontsize',10);
[25836]187 label=getfieldvalue(options,'label','on');
[24001]188
[25836]189 if exist(options,'epsg'),
190 if exist(options,'proj'),
191 error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
192 end
193 proj=epsg2proj(getfieldvalue(options,'epsg'));
194 elseif exist(options,'proj'),
195 proj=getfieldvalue(options,'proj');
196 else
197 proj=epsg2proj(4326);
198 end
199
[24001]200 %read domain:
201 [path,name,ext]=fileparts(self.shpfilename);
202 if ~strcmpi(ext,'shp'),
203 ext='shp';
204 end
205 domain=shpread([self.shppath '/' name '.' ext]);
206
[25836]207 for i=1:length(domain),
208 try,
209 [lat,long] = gdaltransform(domain(i).x,domain(i).y,self.proj,'EPSG:4326');
210 catch me
211 disp(me.message());
212 self.disp();
[24001]213 end
[25836]214 domain(i).x=long; domain(i).y=lat;
215 end
[24001]216
[25836]217 for i=1:length(domain),
[24001]218
[25836]219 %project on x,y,z reference frame.
220 [x,y,z]=AboveGround(domain(i).x,domain(i).y,radius,aboveground);
221 [xt,yt,zt]=AboveGround(domain(i).x+offset,domain(i).y+offset,radius,300000);
222 hold on;
223 if length(x)==1,
224 p=plot3(x,y,z,'k*');
225 set(p,'MarkerSize',markersize);
226 if strcmpi(label,'on'),
[24001]227 t=text(xt,yt,zt,self.shpfilename,'FontSize',fontsize);
[25836]228 end
229 else
230 p=plot3(x,y,z,'k-');
231 mid=floor(length(xt)/2);
232 if strcmpi(label,'on'),
[24001]233 text(xt(mid),yt(mid),zt(mid),self.shpfilename,'FontSize',fontsize);
234 end
235 end
[25836]236 set(p,'Color',color);
237 set(p,'LineWidth',linewidth);
[24001]238
[25836]239 end
[24001]240
241 end % }}}
242 end
243end
Note: See TracBrowser for help on using the repository browser.