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

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

@@ -15 +15 @@
         error('buildgridded error message: gridded not supported for 3d meshes, project on a layer');
 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);
-
-%limits in lat/long
-ullat = max(latlist); ullon = min(lonlist);
-lrlat = min(latlist); lrlon = max(lonlist);
-
-%Find optimal zoom
-if exist(options,'zoom'),
-        zoom = getfieldvalue(options,'zoom');
+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
-        zoom    = optimalzoom(ullat,ullon,lrlat,lrlon);
-        display(['Info: default zoom level ' num2str(zoom)]);
-end
-scale   = 1;
-maxsize = 640;
-transparency = getfieldvalue(options,'transparency',.3);
+        disp('Extracting image from Google maps...');
 
-%convert all these coordinates to pixels
-[ulx, uly]= latlontopixels(ullat, ullon, zoom);
-[lrx, lry]= latlontopixels(lrlat, lrlon, zoom);
+        %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);
 
-%calculate total pixel dimensions of final image
-dx = lrx - ulx;
-dy = uly - lry;
+        %Image corners in lat/long
+        ullat = max(latlist); ullon = min(lonlist);
+        lrlat = min(latlist); lrlon = max(lonlist);
 
-%calculate rows and columns
-cols = ceil(dx/maxsize);
-rows = ceil(dy/maxsize);
+        md=googlemaps(md,ullat,ullon,lrlat,lrlon,options);
+end
 
-%calculate pixel dimensions of each small image
-bottom = 120;
-largura = ceil(dx/cols);
-altura  = ceil(dy/rows);
-alturaplus = altura + bottom;
+%Retrieve image from md
+X = md.radaroverlay.x;
+Y = md.radaroverlay.y;
+final = md.radaroverlay.pwr;
 
-%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
+%Get some options
+transparency = getfieldvalue(options,'transparency',.3);
 
-%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'),
         caxis_opt=getfieldvalue(options,'caxis');
@@ -101 +57 @@
         data_min=min(data_grid(:));
         data_max=max(data_grid(:));
 end
-
 colorm = getcolormap(options);
 image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm);
 
@@ -131 +86 @@
 options=addfielddefault(options,'axis','xy equal off'); % default axis
 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

../trunk-jpl/src/m/plot/googlemaps.m

@@ +1 @@
+function md = googlemaps(md,ullat,ullon,lrlat,lrlon,varargin)
+%GOOGLEMAPS - Extract image from Google maps for given region
+%
+%   Usage:
+%       md = googlemaps(md,ullat,ullon,lrlat,lrlon)
+%       md = googlemaps(md,ullat,ullon,lrlat,lrlon,options)
+%
+%   - ullat,ullon: Upper Left corner latitude and longitude
+%   - lrlat,lrlon: Lower Right corner latitude and longitude
+%
+%   Available options:
+%      - zoom: zoom level, between 1 and 21 (default dynamically calculated)
+
+%Parse inputs
+if nargin<5,
+        help googlemaps
+        error('Wrong usage');
+elseif nargin==5,
+        options=pairoptions;
+else
+        options=varargin{:};
+        if ~isa(options,'pairoptions'),
+                options=pairoptions(varargin{:});
+        end
+end
+
+%Find optimal zoom
+if exist(options,'zoom'),
+        zoom = getfieldvalue(options,'zoom');
+else
+        zoom    = optimalzoom(ullat,ullon,lrlat,lrlon);
+        display(['googlemaps info: default zoom level ' num2str(zoom)]);
+end
+scale   = 1;
+maxsize = 640;
+
+%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 coordinates grids
+[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);
+
+md.radaroverlay.pwr=final;
+md.radaroverlay.x=X;
+md.radaroverlay.y=Y;
+
+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,my]=ll2mercator(lat,lon);
+        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 lon] = mercator2ll(mx,my);
+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