Changeset 15148

Timestamp:

05/29/13 15:50:31 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: separated plot_googlemaps.m and googlemaps.m

Location:

issm/trunk-jpl/src/m/plot

Files:

• googlemaps.m (added)
• plot_googlemaps.m (modified) (3 diffs)

issm/trunk-jpl/src/m/plot/plot_googlemaps.m

@@ -16 +16 @@
 end
 
-%Get xlim and ylim (used to extract radar image)
-xlim=getfieldvalue(options,'xlim',[min(x) max(x)]);
-ylim=getfieldvalue(options,'ylim',[min(y) max(y)]);
-[latlist lonlist]= xy2ll(...
-        [linspace(xlim(1),xlim(2),100) linspace(xlim(2),xlim(2),100) linspace(xlim(2),xlim(1),100) linspace(xlim(1),xlim(1),100)],...
-        [linspace(ylim(1),ylim(1),100) linspace(ylim(1),ylim(2),100) linspace(ylim(2),ylim(2),100) linspace(ylim(2),ylim(1),100)],...
-        +1,45,70);
+if ~any(isnan(md.radaroverlay.x)) & ~any(isnan(md.radaroverlay.y)) & ~any(isnan(md.radaroverlay.pwr)) &...
+                size(md.radaroverlay.pwr,3)==3 & all(size(md.radaroverlay.x)==size(md.radaroverlay.pwr)),
+        disp('plot_googlemaps info: the RGB image held by the model is being used');
+else
+        disp('Extracting image from Google maps...');
 
-%limits in lat/long
-ullat = max(latlist); ullon = min(lonlist);
-lrlat = min(latlist); lrlon = max(lonlist);
+        %Get xlim and ylim (used to extract radar image)
+        xlim=getfieldvalue(options,'xlim',[min(x) max(x)]);
+        ylim=getfieldvalue(options,'ylim',[min(y) max(y)]);
+        [latlist lonlist]= xy2ll(...
+                [linspace(xlim(1),xlim(2),100) linspace(xlim(2),xlim(2),100) linspace(xlim(2),xlim(1),100) linspace(xlim(1),xlim(1),100)],...
+                [linspace(ylim(1),ylim(1),100) linspace(ylim(1),ylim(2),100) linspace(ylim(2),ylim(2),100) linspace(ylim(2),ylim(1),100)],...
+                +1,45,70);
 
-%Find optimal zoom
-if exist(options,'zoom'),
-        zoom = getfieldvalue(options,'zoom');
-else
-        zoom    = optimalzoom(ullat,ullon,lrlat,lrlon);
-        display(['Info: default zoom level ' num2str(zoom)]);
+        %Image corners in lat/long
+        ullat = max(latlist); ullon = min(lonlist);
+        lrlat = min(latlist); lrlon = max(lonlist);
+
+        md=googlemaps(md,ullat,ullon,lrlat,lrlon,options);
 end
-scale   = 1;
-maxsize = 640;
+
+%Retrieve image from md
+X = md.radaroverlay.x;
+Y = md.radaroverlay.y;
+final = md.radaroverlay.pwr;
+
+%Get some options
 transparency = getfieldvalue(options,'transparency',.3);
 
-%convert all these coordinates to pixels
-[ulx, uly]= latlontopixels(ullat, ullon, zoom);
-[lrx, lry]= latlontopixels(lrlat, lrlon, zoom);
-
-%calculate total pixel dimensions of final image
-dx = lrx - ulx;
-dy = uly - lry;
-
-%calculate rows and columns
-cols = ceil(dx/maxsize);
-rows = ceil(dy/maxsize);
-
-%calculate pixel dimensions of each small image
-bottom = 120;
-largura = ceil(dx/cols);
-altura  = ceil(dy/rows);
-alturaplus = altura + bottom;
-
-%Initialize final image
-final = zeros(floor(dy),floor(dx),3);%RGB image
-for x=0:cols-1,
-        for y=0:rows-1,
-                dxn = largura * (0.5 + x);
-                dyn = altura * (0.5 + y);
-                [latn, lonn] = pixelstolatlon(ulx + dxn, uly - dyn - bottom/2, zoom);
-                position = [num2str(latn) ',' num2str(lonn)];
-                disp(['Google Earth tile: ' num2str(x) '/' num2str(cols-1) ' ' num2str(y) '/' num2str(rows-1) ' (center: ' position ')']);
-                params = [...
-                        'center=' position ...
-                        '&zoom=' num2str(zoom)...
-                        '&size=' num2str(largura) 'x' num2str(alturaplus)...
-                        '&maptype=satellite'...
-                        '&sensor=false'...
-                        '&scale=' num2str(scale)];
-                url = ['http://maps.google.com/maps/api/staticmap?' params];
-                [X, map]=imread(url,'png');
-                X=ind2rgb(X,map);
-                indx1 = floor(x*largura)+1;
-                indx2 = min(floor(dx),floor(x*largura)+size(X,2));
-                indy1 = floor(y*altura)+1;
-                indy2 = min(floor(dy),floor(y*altura)+size(X,1));
-                final(indy1:indy2,indx1:indx2,:)=X(1:indy2-indy1+1,1:indx2-indx1+1,:);
-        end
-end
-
-%Create model image
-[gX gY]=meshgrid(ulx:ulx+size(final,2)-1,uly:-1:uly-size(final,1)+1);
-[LAT LON]=pixelstolatlon(gX,gY, zoom);
-[X Y]=ll2xy(LAT,LON,+1,45,70);
+%Prepare grid
 data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X));
-
-%Process data_grid: add white in NaN and correct caxis accordingly
 data_nan=isnan(data_grid);
 if exist(options,'caxis'),
@@ -102 +58 @@
         data_max=max(data_grid(:));
 end
-
 colorm = getcolormap(options);
 image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm);
@@ -132 +87 @@
 applyoptions(md,data,options);
 end
-
-function [px py]=latlontopixels(lat, lon, zoom),
-
-        EARTH_RADIUS = 6378137;
-        EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
-        INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
-        ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
-
-    mx = (lon * ORIGIN_SHIFT) / 180.0;
-    my = log(tan((90 + lat) * pi/360.0))/(pi/180.0);
-    my = (my * ORIGIN_SHIFT) /180.0;
-    res = INITIAL_RESOLUTION / (2^zoom);
-    px = (mx + ORIGIN_SHIFT) / res;
-    py = (my + ORIGIN_SHIFT) / res;
- end
-
-function [lat lon]=pixelstolatlon(px, py, zoom),
-
-        EARTH_RADIUS = 6378137;
-        EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
-        INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
-        ORIGIN_SHIFT = EQUATOR_CIRCUMFERENCE / 2.0;
-
-    res = INITIAL_RESOLUTION / (2^zoom);
-    mx = px * res - ORIGIN_SHIFT;
-    my = py * res - ORIGIN_SHIFT;
-    lat = (my / ORIGIN_SHIFT) * 180.0;
-    lat = 180 / pi * (2*atan(exp(lat*pi/180.0)) - pi/2.0);
-    lon = (mx / ORIGIN_SHIFT) * 180.0;
- end
-function  zoom = optimalzoom(ullat,ullon,lrlat,lrlon)
-
-        EARTH_RADIUS = 6378137;
-        EQUATOR_CIRCUMFERENCE = 2 * pi * EARTH_RADIUS;
-        INITIAL_RESOLUTION = EQUATOR_CIRCUMFERENCE / 256.0;
-
-        optimalsize = 1000; %Number of pixels in final image
-
-        [ulmx ulmy]=ll2mercator(ullat,ullon);
-        [lrmx lrmy]=ll2mercator(lrlat,lrlon);
-        distance = sqrt((lrmx-ulmx)^2 + (lrmy-ulmy)^2);
-
-        zoom1 = floor(log(INITIAL_RESOLUTION*optimalsize/(lrmx-ulmx))/log(2));
-        zoom2 = floor(log(INITIAL_RESOLUTION*optimalsize/(ulmy-lrmy))/log(2));
-
-        zoom=max(zoom1,zoom2);
-
-        zoom = min(max(1,zoom),21);
-end