source: issm/oecreview/Archive/14312-15392/ISSM-15147-15148.diff@ 15393

Last change on this file since 15393 was 15393, checked in by Mathieu Morlighem, 12 years ago

NEW: adding Archive/14312-15392 for oecreview
File size: 10.2 KB
RevLine 
[15393]1Index: ../trunk-jpl/src/m/plot/plot_googlemaps.m
2===================================================================
3--- ../trunk-jpl/src/m/plot/plot_googlemaps.m (revision 15147)
4+++ ../trunk-jpl/src/m/plot/plot_googlemaps.m (revision 15148)
5@@ -15,81 +15,37 @@
6 error('buildgridded error message: gridded not supported for 3d meshes, project on a layer');
7 end
8
9-%Get xlim and ylim (used to extract radar image)
10-xlim=getfieldvalue(options,'xlim',[min(x) max(x)]);
11-ylim=getfieldvalue(options,'ylim',[min(y) max(y)]);
12-[latlist lonlist]= xy2ll(...
13- [linspace(xlim(1),xlim(2),100) linspace(xlim(2),xlim(2),100) linspace(xlim(2),xlim(1),100) linspace(xlim(1),xlim(1),100)],...
14- [linspace(ylim(1),ylim(1),100) linspace(ylim(1),ylim(2),100) linspace(ylim(2),ylim(2),100) linspace(ylim(2),ylim(1),100)],...
15- +1,45,70);
16-
17-%limits in lat/long
18-ullat = max(latlist); ullon = min(lonlist);
19-lrlat = min(latlist); lrlon = max(lonlist);
20-
21-%Find optimal zoom
22-if exist(options,'zoom'),
23- zoom = getfieldvalue(options,'zoom');
24+if ~any(isnan(md.radaroverlay.x)) & ~any(isnan(md.radaroverlay.y)) & ~any(isnan(md.radaroverlay.pwr)) &...
25+ size(md.radaroverlay.pwr,3)==3 & all(size(md.radaroverlay.x)==size(md.radaroverlay.pwr)),
26+ disp('plot_googlemaps info: the RGB image held by the model is being used');
27 else
28- zoom = optimalzoom(ullat,ullon,lrlat,lrlon);
29- display(['Info: default zoom level ' num2str(zoom)]);
30-end
31-scale = 1;
32-maxsize = 640;
33-transparency = getfieldvalue(options,'transparency',.3);
34+ disp('Extracting image from Google maps...');
35
36-%convert all these coordinates to pixels
37-[ulx, uly]= latlontopixels(ullat, ullon, zoom);
38-[lrx, lry]= latlontopixels(lrlat, lrlon, zoom);
39+ %Get xlim and ylim (used to extract radar image)
40+ xlim=getfieldvalue(options,'xlim',[min(x) max(x)]);
41+ ylim=getfieldvalue(options,'ylim',[min(y) max(y)]);
42+ [latlist lonlist]= xy2ll(...
43+ [linspace(xlim(1),xlim(2),100) linspace(xlim(2),xlim(2),100) linspace(xlim(2),xlim(1),100) linspace(xlim(1),xlim(1),100)],...
44+ [linspace(ylim(1),ylim(1),100) linspace(ylim(1),ylim(2),100) linspace(ylim(2),ylim(2),100) linspace(ylim(2),ylim(1),100)],...
45+ +1,45,70);
46
47-%calculate total pixel dimensions of final image
48-dx = lrx - ulx;
49-dy = uly - lry;
50+ %Image corners in lat/long
51+ ullat = max(latlist); ullon = min(lonlist);
52+ lrlat = min(latlist); lrlon = max(lonlist);
53
54-%calculate rows and columns
55-cols = ceil(dx/maxsize);
56-rows = ceil(dy/maxsize);
57+ md=googlemaps(md,ullat,ullon,lrlat,lrlon,options);
58+end
59
60-%calculate pixel dimensions of each small image
61-bottom = 120;
62-largura = ceil(dx/cols);
63-altura = ceil(dy/rows);
64-alturaplus = altura + bottom;
65+%Retrieve image from md
66+X = md.radaroverlay.x;
67+Y = md.radaroverlay.y;
68+final = md.radaroverlay.pwr;
69
70-%Initialize final image
71-final = zeros(floor(dy),floor(dx),3);%RGB image
72-for x=0:cols-1,
73- for y=0:rows-1,
74- dxn = largura * (0.5 + x);
75- dyn = altura * (0.5 + y);
76- [latn, lonn] = pixelstolatlon(ulx + dxn, uly - dyn - bottom/2, zoom);
77- position = [num2str(latn) ',' num2str(lonn)];
78- disp(['Google Earth tile: ' num2str(x) '/' num2str(cols-1) ' ' num2str(y) '/' num2str(rows-1) ' (center: ' position ')']);
79- params = [...
80- 'center=' position ...
81- '&zoom=' num2str(zoom)...
82- '&size=' num2str(largura) 'x' num2str(alturaplus)...
83- '&maptype=satellite'...
84- '&sensor=false'...
85- '&scale=' num2str(scale)];
86- url = ['http://maps.google.com/maps/api/staticmap?' params];
87- [X, map]=imread(url,'png');
88- X=ind2rgb(X,map);
89- indx1 = floor(x*largura)+1;
90- indx2 = min(floor(dx),floor(x*largura)+size(X,2));
91- indy1 = floor(y*altura)+1;
92- indy2 = min(floor(dy),floor(y*altura)+size(X,1));
93- final(indy1:indy2,indx1:indx2,:)=X(1:indy2-indy1+1,1:indx2-indx1+1,:);
94- end
95-end
96+%Get some options
97+transparency = getfieldvalue(options,'transparency',.3);
98
99-%Create model image
100-[gX gY]=meshgrid(ulx:ulx+size(final,2)-1,uly:-1:uly-size(final,1)+1);
101-[LAT LON]=pixelstolatlon(gX,gY, zoom);
102-[X Y]=ll2xy(LAT,LON,+1,45,70);
103+%Prepare grid
104 data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X));
105-
106-%Process data_grid: add white in NaN and correct caxis accordingly
107 data_nan=isnan(data_grid);
108 if exist(options,'caxis'),
109 caxis_opt=getfieldvalue(options,'caxis');
110@@ -101,7 +57,6 @@
111 data_min=min(data_grid(:));
112 data_max=max(data_grid(:));
113 end
114-
115 colorm = getcolormap(options);
116 image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm);
117
118@@ -131,52 +86,3 @@
119 options=addfielddefault(options,'axis','xy equal off'); % default axis
120 applyoptions(md,data,options);
121 end
122-
123-function [px py]=latlontopixels(lat, lon, zoom),
124-
125- EARTH_RADIUS = 6378137;
126- EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
127- INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
128- ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
129-
130- mx = (lon * ORIGIN_SHIFT) / 180.0;
131- my = log(tan((90 + lat) * pi/360.0))/(pi/180.0);
132- my = (my * ORIGIN_SHIFT) /180.0;
133- res = INITIAL_RESOLUTION / (2^zoom);
134- px = (mx + ORIGIN_SHIFT) / res;
135- py = (my + ORIGIN_SHIFT) / res;
136- end
137-
138-function [lat lon]=pixelstolatlon(px, py, zoom),
139-
140- EARTH_RADIUS = 6378137;
141- EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
142- INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
143- ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
144-
145- res = INITIAL_RESOLUTION / (2^zoom);
146- mx = px * res - ORIGIN_SHIFT;
147- my = py * res - ORIGIN_SHIFT;
148- lat = (my / ORIGIN_SHIFT) * 180.0;
149- lat = 180 / pi * (2*atan(exp(lat*pi/180.0)) - pi/2.0);
150- lon = (mx / ORIGIN_SHIFT) * 180.0;
151- end
152-function zoom = optimalzoom(ullat,ullon,lrlat,lrlon)
153-
154- EARTH_RADIUS = 6378137;
155- EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
156- INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
157-
158- optimalsize = 1000; %Number of pixels in final image
159-
160- [ulmx ulmy]=ll2mercator(ullat,ullon);
161- [lrmx lrmy]=ll2mercator(lrlat,lrlon);
162- distance = sqrt((lrmx-ulmx)^2 + (lrmy-ulmy)^2);
163-
164- zoom1 = floor(log(INITIAL_RESOLUTION*optimalsize/(lrmx-ulmx))/log(2));
165- zoom2 = floor(log(INITIAL_RESOLUTION*optimalsize/(ulmy-lrmy))/log(2));
166-
167- zoom=max(zoom1,zoom2);
168-
169- zoom = min(max(1,zoom),21);
170-end
171Index: ../trunk-jpl/src/m/plot/googlemaps.m
172===================================================================
173--- ../trunk-jpl/src/m/plot/googlemaps.m (revision 0)
174+++ ../trunk-jpl/src/m/plot/googlemaps.m (revision 15148)
175@@ -0,0 +1,133 @@
176+function md = googlemaps(md,ullat,ullon,lrlat,lrlon,varargin)
177+%GOOGLEMAPS - Extract image from Google maps for given region
178+%
179+% Usage:
180+% md = googlemaps(md,ullat,ullon,lrlat,lrlon)
181+% md = googlemaps(md,ullat,ullon,lrlat,lrlon,options)
182+%
183+% - ullat,ullon: Upper Left corner latitude and longitude
184+% - lrlat,lrlon: Lower Right corner latitude and longitude
185+%
186+% Available options:
187+% - zoom: zoom level, between 1 and 21 (default dynamically calculated)
188+
189+%Parse inputs
190+if nargin<5,
191+ help googlemaps
192+ error('Wrong usage');
193+elseif nargin==5,
194+ options=pairoptions;
195+else
196+ options=varargin{:};
197+ if ~isa(options,'pairoptions'),
198+ options=pairoptions(varargin{:});
199+ end
200+end
201+
202+%Find optimal zoom
203+if exist(options,'zoom'),
204+ zoom = getfieldvalue(options,'zoom');
205+else
206+ zoom = optimalzoom(ullat,ullon,lrlat,lrlon);
207+ display(['googlemaps info: default zoom level ' num2str(zoom)]);
208+end
209+scale = 1;
210+maxsize = 640;
211+
212+%convert all these coordinates to pixels
213+[ulx, uly]= latlontopixels(ullat, ullon, zoom);
214+[lrx, lry]= latlontopixels(lrlat, lrlon, zoom);
215+
216+%calculate total pixel dimensions of final image
217+dx = lrx - ulx;
218+dy = uly - lry;
219+
220+%calculate rows and columns
221+cols = ceil(dx/maxsize);
222+rows = ceil(dy/maxsize);
223+
224+%calculate pixel dimensions of each small image
225+bottom = 120;
226+largura = ceil(dx/cols);
227+altura = ceil(dy/rows);
228+alturaplus = altura + bottom;
229+
230+%Initialize final image
231+final = zeros(floor(dy),floor(dx),3);%RGB image
232+for x=0:cols-1,
233+ for y=0:rows-1,
234+ dxn = largura * (0.5 + x);
235+ dyn = altura * (0.5 + y);
236+ [latn, lonn] = pixelstolatlon(ulx + dxn, uly - dyn - bottom/2, zoom);
237+ position = [num2str(latn) ',' num2str(lonn)];
238+ disp(['Google Earth tile: ' num2str(x) '/' num2str(cols-1) ' ' num2str(y) '/' num2str(rows-1) ' (center: ' position ')']);
239+ params = [...
240+ 'center=' position ...
241+ '&zoom=' num2str(zoom)...
242+ '&size=' num2str(largura) 'x' num2str(alturaplus)...
243+ '&maptype=satellite'...
244+ '&sensor=false'...
245+ '&scale=' num2str(scale)];
246+ url = ['http://maps.google.com/maps/api/staticmap?' params];
247+ [X, map]=imread(url,'png');
248+ X=ind2rgb(X,map);
249+ indx1 = floor(x*largura)+1;
250+ indx2 = min(floor(dx),floor(x*largura)+size(X,2));
251+ indy1 = floor(y*altura)+1;
252+ indy2 = min(floor(dy),floor(y*altura)+size(X,1));
253+ final(indy1:indy2,indx1:indx2,:)=X(1:indy2-indy1+1,1:indx2-indx1+1,:);
254+ end
255+end
256+
257+%Create coordinates grids
258+[gX gY]=meshgrid(ulx:ulx+size(final,2)-1,uly:-1:uly-size(final,1)+1);
259+[LAT LON]=pixelstolatlon(gX,gY, zoom);
260+[X Y]=ll2xy(LAT,LON,+1,45,70);
261+
262+md.radaroverlay.pwr=final;
263+md.radaroverlay.x=X;
264+md.radaroverlay.y=Y;
265+
266+end
267+function [px py]=latlontopixels(lat, lon, zoom),
268+ EARTH_RADIUS = 6378137;
269+ EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
270+ INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
271+ ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
272+
273+ [mx,my]=ll2mercator(lat,lon);
274+ res = INITIAL_RESOLUTION / (2^zoom);
275+ px = (mx + ORIGIN_SHIFT) / res;
276+ py = (my + ORIGIN_SHIFT) / res;
277+end
278+
279+function [lat lon]=pixelstolatlon(px, py, zoom),
280+ EARTH_RADIUS = 6378137;
281+ EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
282+ INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
283+ ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
284+
285+ res = INITIAL_RESOLUTION / (2^zoom);
286+ mx = px * res - ORIGIN_SHIFT;
287+ my = py * res - ORIGIN_SHIFT;
288+ [lat lon] = mercator2ll(mx,my);
289+end
290+function zoom = optimalzoom(ullat,ullon,lrlat,lrlon)
291+
292+ EARTH_RADIUS = 6378137;
293+ EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
294+ INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
295+
296+ optimalsize = 1000; %Number of pixels in final image
297+
298+ [ulmx ulmy]=ll2mercator(ullat,ullon);
299+ [lrmx lrmy]=ll2mercator(lrlat,lrlon);
300+ distance = sqrt((lrmx-ulmx)^2 + (lrmy-ulmy)^2);
301+
302+ zoom1 = floor(log(INITIAL_RESOLUTION*optimalsize/(lrmx-ulmx))/log(2));
303+ zoom2 = floor(log(INITIAL_RESOLUTION*optimalsize/(ulmy-lrmy))/log(2));
304+
305+ zoom=max(zoom1,zoom2);
306+
307+ zoom = min(max(1,zoom),21);
308+end
Note: See TracBrowser for help on using the repository browser.