Index: /issm/trunk-jpl/src/m/plot/googlemaps.m =================================================================== --- /issm/trunk-jpl/src/m/plot/googlemaps.m (revision 15148) +++ /issm/trunk-jpl/src/m/plot/googlemaps.m (revision 15148) @@ -0,0 +1,133 @@ +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 Index: /issm/trunk-jpl/src/m/plot/plot_googlemaps.m =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_googlemaps.m (revision 15147) +++ /issm/trunk-jpl/src/m/plot/plot_googlemaps.m (revision 15148) @@ -16,79 +16,35 @@ 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,5 +58,4 @@ 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,51 +87,2 @@ 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