source: issm/trunk-jpl/src/py3/plot/plot_googlemaps.py@ 23677

function plot_googlemaps(md,data,options,plotlines,plotcols,i)
%PLOT_GOOGLEMAPS - superimpose Google maps to a given field
%
%   Usage:
%      plot_googlemaps(md,data,options,plotlines,plotcols,i)
%
%   See also: PLOTMODEL

%process mesh and data
[x y z elements is2d isplanet]=processmesh(md,[],options);
[data datatype]=processdata(md,data,options);

%check is2d
if ~is2d, 
        error('buildgridded error message: gridded not supported for 3d meshes, project on a layer');
end

if ~any(isnan(md.radaroverlay.x(:))) & ~any(isnan(md.radaroverlay.y(:))) & ~any(isnan(md.radaroverlay.pwr(:))) ...
                & size(md.radaroverlay.pwr,3)==3 & size(md.radaroverlay.x,2)==size(md.radaroverlay.pwr,2),
        disp('plot_googlemaps info: the RGB image held by the model is being used');
else
        disp('Extracting image from Google maps...');

        %Get xlim and ylim (used to extract radar image)
        xlim=getfieldvalue(options,'xlim',[min(x) max(x)])/getfieldvalue(options,'unit',1);
        ylim=getfieldvalue(options,'ylim',[min(y) max(y)])/getfieldvalue(options,'unit',1);
        if md.mesh.epsg==3413, %UPS Greenland
                [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);
        elseif md.mesh.epsg==3031, %UPS Antarctica
                [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,0,71);
        elseif md.mesh.epsg==26906, %UTM 6V Columbia Glacier Alaska
                [latlist lonlist]= utm2ll(...
                        [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)],...
                        6);
        elseif numel(md.mesh.lat)==numel(md.mesh.x),
                latlist = md.mesh.lat; %That might work?
                lonlist = md.mesh.long;
        else
                error('EPSG code not supported yet, and no lat long found in md.mesh');
        end

        %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

%Process image from model
final = double(md.radaroverlay.pwr)/double(max(md.radaroverlay.pwr(:))); %rescale between 0 and 1

%Get some options
transparency = getfieldvalue(options,'transparency',.3);

%Prepare grid
if size(md.radaroverlay.x,1)==1 | size(md.radaroverlay.x,2)==1,
        x_m = md.radaroverlay.x;
        y_m = md.radaroverlay.y;
        data_grid=InterpFromMeshToGrid(elements,x/getfieldvalue(options,'unit',1),y/getfieldvalue(options,'unit',1),data,x_m,y_m,NaN);
else
        X = md.radaroverlay.x;
        Y = md.radaroverlay.y;
        data_grid=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X(:),Y(:),'default',NaN); data_grid=reshape(data_grid,size(X));
        x_m=X(1,:); y_m=Y(:,1);
end
data_nan=isnan(data_grid);
if exist(options,'caxis'),
        caxis_opt=getfieldvalue(options,'caxis');
        data_grid(find(data_grid<caxis_opt(1)))=caxis_opt(1);
        data_grid(find(data_grid>caxis_opt(2)))=caxis_opt(2);
        data_min=caxis_opt(1);
        data_max=caxis_opt(2);
else
        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);

if exist(options,'shaded'),
        a    = -45;
        scut = 0.2;
        c    = 1;
        % computes lighting from elevation gradient
        [fx,fy] = gradient(data_grid,x_m,y_m);
        fxy = -fx*sind(a) - fy*cosd(a);
        clear fx fy % free some memory...
        fxy(isnan(fxy)) = 0;

        % computes maximum absolute gradient (median-style), normalizes, saturates and duplicates in 3-D matrix
        r = repmat(max(min(fxy/nmedian(abs(fxy),1 - scut/100),1),-1),[1,1,3]);

        % applies contrast using exponent
        rp = (1 - abs(r)).^c;
        image_rgb = image_rgb.*rp;

        % lighter for positive gradient
        k = find(r > 0);
        image_rgb(k) = image_rgb(k) + (1 - rp(k));
end

alpha=ones(size(data_grid));
alpha(find(~data_nan))=transparency;
alpha=repmat(alpha,[1 1 3]);

final=alpha.*final+(1-alpha).*image_rgb;

%Select plot area 
subplotmodel(plotlines,plotcols,i,options);

h=imagesc(x_m*getfieldvalue(options,'unit',1),y_m*getfieldvalue(options,'unit',1),final);

%last step: mesh gridded?
if exist(options,'edgecolor'),
        A=elements(:,1); B=elements(:,2); C=elements(:,3); 
        patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor'));
end

%Apply options
if ~isnan(data_min),
        options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready
end
options=addfielddefault(options,'axis','xy equal off'); % default axis
applyoptions(md,data,options);
end