02.03.2016 Google Maps в MATLAB
Материал из SRNS
Добавить Google-карту в свою MATLAB-модель? Легко!
В community Mathworks нашлась функция plot_google_map:
plot_google_map
function varargout = plot_google_map(varargin)
% function h = plot_google_map(varargin)
% Plots a google map on the current axes using the Google Static Maps API
%
% USAGE:
% h = plot_google_map(Property, Value,...)
% Plots the map on the given axes. Used also if no output is specified
%
% Or:
% [lonVec latVec imag] = plot_google_map(Property, Value,...)
% Returns the map without plotting it
%
% PROPERTIES:
% Axis - Axis handle. If not given, gca is used. (LP)
% Height (640) - Height of the image in pixels (max 640)
% Width (640) - Width of the image in pixels (max 640)
% Scale (2) - (1/2) Resolution scale factor. Using Scale=2 will
% double the resulotion of the downloaded image (up
% to 1280x1280) and will result in finer rendering,
% but processing time will be longer.
% MapType - ('roadmap') Type of map to return. Any of [roadmap,
% satellite, terrain, hybrid]. See the Google Maps API for
% more information.
% Alpha (1) - (0-1) Transparency level of the map (0 is fully
% transparent). While the map is always moved to the
% bottom of the plot (i.e. will not hide previously
% drawn items), this can be useful in order to increase
% readability if many colors are plotted
% (using SCATTER for example).
% ShowLabels (1) - (0/1) Controls whether to display city/street textual labels on the map
% Language - (string) A 2 letter ISO 639-1 language code for displaying labels in a
% local language instead of English (where available).
% For example, for Chinese use:
% plot_google_map('language','zh')
% For the list of codes, see:
% http://en.wikipedia.org/wiki/List_of_ISO_639-1_codes
% Marker - The marker argument is a text string with fields
% conforming to the Google Maps API. The
% following are valid examples:
% '43.0738740,-70.713993' (default midsize orange marker)
% '43.0738740,-70.713993,blue' (midsize blue marker)
% '43.0738740,-70.713993,yellowa' (midsize yellow
% marker with label "A")
% '43.0738740,-70.713993,tinyredb' (tiny red marker
% with label "B")
% Refresh (1) - (0/1) defines whether to automatically refresh the
% map upon zoom/pan action on the figure.
% AutoAxis (1) - (0/1) defines whether to automatically adjust the axis
% of the plot to avoid the map being stretched.
% This will adjust the span to be correct
% according to the shape of the map axes.
% FigureResizeUpdate (1) - (0/1) defines whether to automatically refresh the
% map upon resizing the figure. This will ensure map
% isn't stretched after figure resize.
% APIKey - (string) set your own API key which you obtained from Google:
% http://developers.google.com/maps/documentation/staticmaps/#api_key
% This will enable up to 25,000 map requests per day,
% compared to a few hundred requests without a key.
% To set the key, use:
% plot_google_map('APIKey','SomeLongStringObtaindFromGoogle')
% You need to do this only once to set the key.
% To disable the use of a key, use:
% plot_google_map('APIKey','')
%
% OUTPUT:
% h - Handle to the plotted map
%
% lonVect - Vector of Longidute coordinates (WGS84) of the image
% latVect - Vector of Latidute coordinates (WGS84) of the image
% imag - Image matrix (height,width,3) of the map
%
% EXAMPLE - plot a map showing some capitals in Europe:
% lat = [48.8708 51.5188 41.9260 40.4312 52.523 37.982];
% lon = [2.4131 -0.1300 12.4951 -3.6788 13.415 23.715];
% plot(lon,lat,'.r','MarkerSize',20)
% plot_google_map
%
% References:
% http://www.mathworks.com/matlabcentral/fileexchange/24113
% http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/
% http://developers.google.com/maps/documentation/staticmaps/
%
% Acknowledgements:
% Val Schmidt for his submission of get_google_map.m
%
% Author:
% Zohar Bar-Yehuda
%
% Version 1.6 - 12/11/2015
% - Use system temp folder for writing image files (with fallback to current dir if missing write permissions)
% Version 1.5 - 20/11/2014
% - Support for MATLAB R2014b
% - several fixes complex layouts: several maps in one figure,
% map inside a panel, specifying axis handle as input (thanks to Luke Plausin)
% Version 1.4 - 25/03/2014
% - Added the language parameter for showing labels in a local language
% - Display the URL on error to allow easier debugging of API errors
% Version 1.3 - 06/10/2013
% - Improved functionality of AutoAxis, which now handles any shape of map axes.
% Now also updates the extent of the map if the figure is resized.
% - Added the showLabels parameter which allows hiding the textual labels on the map.
% Version 1.2 - 16/06/2012
% - Support use of the "scale=2" parameter by default for finer rendering (set scale=1 if too slow).
% - Auto-adjust axis extent so the map isn't stretched.
% - Set and use an API key which enables a much higher usage volume per day.
% Version 1.1 - 25/08/2011
persistent apiKey useTemp
if isnumeric(apiKey)
% first run, check if API key file exists
if exist('api_key.mat','file')
load api_key
else
apiKey = '';
end
end
if isempty(useTemp)
% first run, check we we have wrtie access to the temp folder
try
tempfilename = tempname;
fid = fopen(tempfilename, 'w');
if fid > 0
fclose(fid);
useTemp = true;
delete(tempfilename);
else
% Don't have write access to temp folder or it doesn't exist, fallback to current dir
useTemp = false;
end
catch
% in case tempname fails for some reason
useTemp = false;
end
end
hold on
% Default parametrs
axHandle = gca;
height = 640;
width = 640;
scale = 2;
maptype = 'roadmap';
alphaData = 1;
autoRefresh = 1;
figureResizeUpdate = 1;
autoAxis = 1;
showLabels = 1;
language = '';
markeridx = 1;
markerlist = {};
% Handle input arguments
if nargin >= 2
for idx = 1:2:length(varargin)
switch lower(varargin{idx})
case 'axis'
axHandle = varargin{idx+1};
case 'height'
height = varargin{idx+1};
case 'width'
width = varargin{idx+1};
case 'maptype'
maptype = varargin{idx+1};
case 'alpha'
alphaData = varargin{idx+1};
case 'refresh'
autoRefresh = varargin{idx+1};
case 'showlabels'
showLabels = varargin{idx+1};
case 'figureresizeupdate'
figureResizeUpdate = varargin{idx+1};
case 'language'
language = varargin{idx+1};
case 'marker'
markerlist{markeridx} = varargin{idx+1};
markeridx = markeridx + 1;
case 'autoaxis'
autoAxis = varargin{idx+1};
case 'apikey'
apiKey = varargin{idx+1}; % set new key
% save key to file
funcFile = which('plot_google_map.m');
pth = fileparts(funcFile);
keyFile = fullfile(pth,'api_key.mat');
save(keyFile,'apiKey')
otherwise
error(['Unrecognized variable: ' varargin{idx}])
end
end
end
if height > 640
height = 640;
end
if width > 640
width = 640;
end
% Store paramters in axis handle (for auto refresh callbacks)
ud = get(axHandle, 'UserData');
if isempty(ud)
% explicitly set as struct to avoid warnings
ud = struct;
end
ud.gmap_params = varargin;
set(axHandle, 'UserData', ud);
curAxis = axis(axHandle);
% Enforce Latitude constraints of EPSG:900913
if curAxis(3) < -85
curAxis(3) = -85;
end
if curAxis(4) > 85
curAxis(4) = 85;
end
% Enforce longitude constrains
if curAxis(1) < -180
curAxis(1) = -180;
end
if curAxis(1) > 180
curAxis(1) = 0;
end
if curAxis(2) > 180
curAxis(2) = 180;
end
if curAxis(2) < -180
curAxis(2) = 0;
end
if isequal(curAxis,[0 1 0 1]) % probably an empty figure
% display world map
curAxis = [-200 200 -85 85];
axis(curAxis)
end
if autoAxis
% adjust current axis limit to avoid strectched maps
[xExtent,yExtent] = latLonToMeters(curAxis(3:4), curAxis(1:2) );
xExtent = diff(xExtent); % just the size of the span
yExtent = diff(yExtent);
% get axes aspect ratio
drawnow
org_units = get(axHandle,'Units');
set(axHandle,'Units','Pixels')
ax_position = get(axHandle,'position');
set(axHandle,'Units',org_units)
aspect_ratio = ax_position(4) / ax_position(3);
if xExtent*aspect_ratio > yExtent
centerX = mean(curAxis(1:2));
centerY = mean(curAxis(3:4));
spanX = (curAxis(2)-curAxis(1))/2;
spanY = (curAxis(4)-curAxis(3))/2;
% enlarge the Y extent
spanY = spanY*xExtent*aspect_ratio/yExtent; % new span
if spanY > 85
spanX = spanX * 85 / spanY;
spanY = spanY * 85 / spanY;
end
curAxis(1) = centerX-spanX;
curAxis(2) = centerX+spanX;
curAxis(3) = centerY-spanY;
curAxis(4) = centerY+spanY;
elseif yExtent > xExtent*aspect_ratio
centerX = mean(curAxis(1:2));
centerY = mean(curAxis(3:4));
spanX = (curAxis(2)-curAxis(1))/2;
spanY = (curAxis(4)-curAxis(3))/2;
% enlarge the X extent
spanX = spanX*yExtent/(xExtent*aspect_ratio); % new span
if spanX > 180
spanY = spanY * 180 / spanX;
spanX = spanX * 180 / spanX;
end
curAxis(1) = centerX-spanX;
curAxis(2) = centerX+spanX;
curAxis(3) = centerY-spanY;
curAxis(4) = centerY+spanY;
end
% Enforce Latitude constraints of EPSG:900913
if curAxis(3) < -85
curAxis(3:4) = curAxis(3:4) + (-85 - curAxis(3));
end
if curAxis(4) > 85
curAxis(3:4) = curAxis(3:4) + (85 - curAxis(4));
end
axis(axHandle, curAxis); % update axis as quickly as possible, before downloading new image
drawnow
end
% Delete previous map from plot (if exists)
if nargout <= 1 % only if in plotting mode
curChildren = get(axHandle,'children');
map_objs = findobj(curChildren,'tag','gmap');
bd_callback = [];
for idx = 1:length(map_objs)
if ~isempty(get(map_objs(idx),'ButtonDownFcn'))
% copy callback properties from current map
bd_callback = get(map_objs(idx),'ButtonDownFcn');
end
end
delete(map_objs)
end
% Calculate zoom level for current axis limits
[xExtent,yExtent] = latLonToMeters(curAxis(3:4), curAxis(1:2) );
minResX = diff(xExtent) / width;
minResY = diff(yExtent) / height;
minRes = max([minResX minResY]);
tileSize = 256;
initialResolution = 2 * pi * 6378137 / tileSize; % 156543.03392804062 for tileSize 256 pixels
zoomlevel = floor(log2(initialResolution/minRes));
% Enforce valid zoom levels
if zoomlevel < 0
zoomlevel = 0;
end
if zoomlevel > 19
zoomlevel = 19;
end
% Calculate center coordinate in WGS1984
lat = (curAxis(3)+curAxis(4))/2;
lon = (curAxis(1)+curAxis(2))/2;
% Construct query URL
preamble = 'http://maps.googleapis.com/maps/api/staticmap';
location = ['?center=' num2str(lat,10) ',' num2str(lon,10)];
zoomStr = ['&zoom=' num2str(zoomlevel)];
sizeStr = ['&scale=' num2str(scale) '&size=' num2str(width) 'x' num2str(height)];
maptypeStr = ['&maptype=' maptype ];
if ~isempty(apiKey)
keyStr = ['&key=' apiKey];
else
keyStr = '';
end
markers = '&markers=';
for idx = 1:length(markerlist)
if idx < length(markerlist)
markers = [markers markerlist{idx} '%7C'];
else
markers = [markers markerlist{idx}];
end
end
if showLabels == 0
labelsStr = '&style=feature:all|element:labels|visibility:off';
else
labelsStr = '';
end
if ~isempty(language)
languageStr = ['&language=' language];
else
languageStr = '';
end
if ismember(maptype,{'satellite','hybrid'})
filename = 'tmp.jpg';
format = '&format=jpg';
convertNeeded = 0;
else
filename = 'tmp.png';
format = '&format=png';
convertNeeded = 1;
end
sensor = '&sensor=false';
url = [preamble location zoomStr sizeStr maptypeStr format markers labelsStr languageStr sensor keyStr];
% Get the image
if useTemp
filepath = fullfile(tempdir, filename);
else
filepath = filename;
end
try
urlwrite(url,filepath);
catch % error downloading map
warning(['Unable to download map form Google Servers.\n' ...
'Matlab error was: %s\n\n' ...
'Possible reasons: missing write permissions, no network connection, quota exceeded, or some other error.\n' ...
'Consider using an API key if quota problems persist.\n\n' ...
'To debug, try pasting the following URL in your browser, which may result in a more informative error:\n%s'], lasterr, url);
varargout{1} = [];
varargout{2} = [];
varargout{3} = [];
return
end
[M Mcolor] = imread(filepath);
M = cast(M,'double');
delete(filepath); % delete temp file
width = size(M,2);
height = size(M,1);
% Calculate a meshgrid of pixel coordinates in EPSG:900913
centerPixelY = round(height/2);
centerPixelX = round(width/2);
[centerX,centerY] = latLonToMeters(lat, lon ); % center coordinates in EPSG:900913
curResolution = initialResolution / 2^zoomlevel/scale; % meters/pixel (EPSG:900913)
xVec = centerX + ((1:width)-centerPixelX) * curResolution; % x vector
yVec = centerY + ((height:-1:1)-centerPixelY) * curResolution; % y vector
[xMesh,yMesh] = meshgrid(xVec,yVec); % construct meshgrid
% convert meshgrid to WGS1984
[lonMesh,latMesh] = metersToLatLon(xMesh,yMesh);
% We now want to convert the image from a colormap image with an uneven
% mesh grid, into an RGB truecolor image with a uniform grid.
% This would enable displaying it with IMAGE, instead of PCOLOR.
% Advantages are:
% 1) faster rendering
% 2) makes it possible to display together with other colormap annotations (PCOLOR, SCATTER etc.)
% Convert image from colormap type to RGB truecolor (if PNG is used)
if convertNeeded
imag = zeros(height,width,3);
for idx = 1:3
imag(:,:,idx) = reshape(Mcolor(M(:)+1+(idx-1)*size(Mcolor,1)),height,width);
end
else
imag = M/255;
end
% Next, project the data into a uniform WGS1984 grid
sizeFactor = 1; % factoring of new image
uniHeight = round(height*sizeFactor);
uniWidth = round(width*sizeFactor);
latVect = linspace(latMesh(1,1),latMesh(end,1),uniHeight);
lonVect = linspace(lonMesh(1,1),lonMesh(1,end),uniWidth);
[uniLonMesh,uniLatMesh] = meshgrid(lonVect,latVect);
uniImag = zeros(uniHeight,uniWidth,3);
% old version (projection using INTERP2)
% for idx = 1:3
% % 'nearest' method is the fastest. difference from other methods is neglible
% uniImag(:,:,idx) = interp2(lonMesh,latMesh,imag(:,:,idx),uniLonMesh,uniLatMesh,'nearest');
% end
uniImag = myTurboInterp2(lonMesh,latMesh,imag,uniLonMesh,uniLatMesh);
if nargout <= 1 % plot map
% display image
hold(axHandle, 'on');
cax = caxis;
h = image(lonVect,latVect,uniImag, 'Parent', axHandle);
caxis(cax); % Preserve caxis that is sometimes changed by the call to image()
set(axHandle,'YDir','Normal')
set(h,'tag','gmap')
set(h,'AlphaData',alphaData)
% add a dummy image to allow pan/zoom out to x2 of the image extent
h_tmp = image(lonVect([1 end]),latVect([1 end]),zeros(2),'Visible','off', 'Parent', axHandle);
set(h_tmp,'tag','gmap')
% older version (display without conversion to uniform grid)
% h =pcolor(lonMesh,latMesh,(M));
% colormap(Mcolor)
% caxis([0 255])
% warning off % to avoid strange rendering warnings
% shading flat
uistack(h,'bottom') % move map to bottom (so it doesn't hide previously drawn annotations)
axis(axHandle, curAxis) % restore original zoom
if nargout == 1
varargout{1} = h;
end
% if auto-refresh mode - override zoom callback to allow autumatic
% refresh of map upon zoom actions.
figHandle = axHandle;
while ~strcmpi(get(figHandle, 'Type'), 'figure')
% Recursively search for parent figure in case axes are in a panel
figHandle = get(figHandle, 'Parent');
end
zoomHandle = zoom(axHandle);
panHandle = pan(figHandle); % This isn't ideal, doesn't work for contained axis
if autoRefresh
set(zoomHandle,'ActionPostCallback',@update_google_map);
set(panHandle, 'ActionPostCallback', @update_google_map);
else % disable zoom override
set(zoomHandle,'ActionPostCallback',[]);
set(panHandle, 'ActionPostCallback',[]);
end
% set callback for figure resize function, to update extents if figure
% is streched.
if figureResizeUpdate &&isempty(get(figHandle, 'ResizeFcn'))
% set only if not already set by someone else
set(figHandle, 'ResizeFcn', @update_google_map_fig);
end
% set callback properties
set(h,'ButtonDownFcn',bd_callback);
else % don't plot, only return map
varargout{1} = lonVect;
varargout{2} = latVect;
varargout{3} = uniImag;
end
% Coordinate transformation functions
function [lon,lat] = metersToLatLon(x,y)
% Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum
originShift = 2 * pi * 6378137 / 2.0; % 20037508.342789244
lon = (x ./ originShift) * 180;
lat = (y ./ originShift) * 180;
lat = 180 / pi * (2 * atan( exp( lat * pi / 180)) - pi / 2);
function [x,y] = latLonToMeters(lat, lon )
% Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913"
originShift = 2 * pi * 6378137 / 2.0; % 20037508.342789244
x = lon * originShift / 180;
y = log(tan((90 + lat) * pi / 360 )) / (pi / 180);
y = y * originShift / 180;
function ZI = myTurboInterp2(X,Y,Z,XI,YI)
% An extremely fast nearest neighbour 2D interpolation, assuming both input
% and output grids consist only of squares, meaning:
% - uniform X for each column
% - uniform Y for each row
XI = XI(1,:);
X = X(1,:);
YI = YI(:,1);
Y = Y(:,1);
xiPos = nan*ones(size(XI));
xLen = length(X);
yiPos = nan*ones(size(YI));
yLen = length(Y);
% find x conversion
xPos = 1;
for idx = 1:length(xiPos)
if XI(idx) >= X(1) && XI(idx) <= X(end)
while xPos < xLen && X(xPos+1)<XI(idx)
xPos = xPos + 1;
end
diffs = abs(X(xPos:xPos+1)-XI(idx));
if diffs(1) < diffs(2)
xiPos(idx) = xPos;
else
xiPos(idx) = xPos + 1;
end
end
end
% find y conversion
yPos = 1;
for idx = 1:length(yiPos)
if YI(idx) <= Y(1) && YI(idx) >= Y(end)
while yPos < yLen && Y(yPos+1)>YI(idx)
yPos = yPos + 1;
end
diffs = abs(Y(yPos:yPos+1)-YI(idx));
if diffs(1) < diffs(2)
yiPos(idx) = yPos;
else
yiPos(idx) = yPos + 1;
end
end
end
ZI = Z(yiPos,xiPos,:);
function update_google_map(obj,evd)
% callback function for auto-refresh
drawnow;
try
axHandle = evd.Axes;
catch ex
% Event doesn't contain the correct axes. Panic!
axHandle = gca;
end
ud = get(axHandle, 'UserData');
if isfield(ud, 'gmap_params')
params = ud.gmap_params;
plot_google_map(params{:});
end
function update_google_map_fig(obj,evd)
% callback function for auto-refresh
drawnow;
axes_objs = findobj(get(gcf,'children'),'type','axes');
for idx = 1:length(axes_objs)
if ~isempty(findobj(get(axes_objs(idx),'children'),'tag','gmap'));
ud = get(axes_objs(idx), 'UserData');
if isfield(ud, 'gmap_params')
params = ud.gmap_params;
else
params = {};
end
% Add axes to inputs if needed
if ~sum(strcmpi(params, 'Axis'))
params = [params, {'Axis', axes_objs(idx)}];
end
plot_google_map(params{:});
end
end
% function h = plot_google_map(varargin)
% Plots a google map on the current axes using the Google Static Maps API
%
% USAGE:
% h = plot_google_map(Property, Value,...)
% Plots the map on the given axes. Used also if no output is specified
%
% Or:
% [lonVec latVec imag] = plot_google_map(Property, Value,...)
% Returns the map without plotting it
%
% PROPERTIES:
% Axis - Axis handle. If not given, gca is used. (LP)
% Height (640) - Height of the image in pixels (max 640)
% Width (640) - Width of the image in pixels (max 640)
% Scale (2) - (1/2) Resolution scale factor. Using Scale=2 will
% double the resulotion of the downloaded image (up
% to 1280x1280) and will result in finer rendering,
% but processing time will be longer.
% MapType - ('roadmap') Type of map to return. Any of [roadmap,
% satellite, terrain, hybrid]. See the Google Maps API for
% more information.
% Alpha (1) - (0-1) Transparency level of the map (0 is fully
% transparent). While the map is always moved to the
% bottom of the plot (i.e. will not hide previously
% drawn items), this can be useful in order to increase
% readability if many colors are plotted
% (using SCATTER for example).
% ShowLabels (1) - (0/1) Controls whether to display city/street textual labels on the map
% Language - (string) A 2 letter ISO 639-1 language code for displaying labels in a
% local language instead of English (where available).
% For example, for Chinese use:
% plot_google_map('language','zh')
% For the list of codes, see:
% http://en.wikipedia.org/wiki/List_of_ISO_639-1_codes
% Marker - The marker argument is a text string with fields
% conforming to the Google Maps API. The
% following are valid examples:
% '43.0738740,-70.713993' (default midsize orange marker)
% '43.0738740,-70.713993,blue' (midsize blue marker)
% '43.0738740,-70.713993,yellowa' (midsize yellow
% marker with label "A")
% '43.0738740,-70.713993,tinyredb' (tiny red marker
% with label "B")
% Refresh (1) - (0/1) defines whether to automatically refresh the
% map upon zoom/pan action on the figure.
% AutoAxis (1) - (0/1) defines whether to automatically adjust the axis
% of the plot to avoid the map being stretched.
% This will adjust the span to be correct
% according to the shape of the map axes.
% FigureResizeUpdate (1) - (0/1) defines whether to automatically refresh the
% map upon resizing the figure. This will ensure map
% isn't stretched after figure resize.
% APIKey - (string) set your own API key which you obtained from Google:
% http://developers.google.com/maps/documentation/staticmaps/#api_key
% This will enable up to 25,000 map requests per day,
% compared to a few hundred requests without a key.
% To set the key, use:
% plot_google_map('APIKey','SomeLongStringObtaindFromGoogle')
% You need to do this only once to set the key.
% To disable the use of a key, use:
% plot_google_map('APIKey','')
%
% OUTPUT:
% h - Handle to the plotted map
%
% lonVect - Vector of Longidute coordinates (WGS84) of the image
% latVect - Vector of Latidute coordinates (WGS84) of the image
% imag - Image matrix (height,width,3) of the map
%
% EXAMPLE - plot a map showing some capitals in Europe:
% lat = [48.8708 51.5188 41.9260 40.4312 52.523 37.982];
% lon = [2.4131 -0.1300 12.4951 -3.6788 13.415 23.715];
% plot(lon,lat,'.r','MarkerSize',20)
% plot_google_map
%
% References:
% http://www.mathworks.com/matlabcentral/fileexchange/24113
% http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/
% http://developers.google.com/maps/documentation/staticmaps/
%
% Acknowledgements:
% Val Schmidt for his submission of get_google_map.m
%
% Author:
% Zohar Bar-Yehuda
%
% Version 1.6 - 12/11/2015
% - Use system temp folder for writing image files (with fallback to current dir if missing write permissions)
% Version 1.5 - 20/11/2014
% - Support for MATLAB R2014b
% - several fixes complex layouts: several maps in one figure,
% map inside a panel, specifying axis handle as input (thanks to Luke Plausin)
% Version 1.4 - 25/03/2014
% - Added the language parameter for showing labels in a local language
% - Display the URL on error to allow easier debugging of API errors
% Version 1.3 - 06/10/2013
% - Improved functionality of AutoAxis, which now handles any shape of map axes.
% Now also updates the extent of the map if the figure is resized.
% - Added the showLabels parameter which allows hiding the textual labels on the map.
% Version 1.2 - 16/06/2012
% - Support use of the "scale=2" parameter by default for finer rendering (set scale=1 if too slow).
% - Auto-adjust axis extent so the map isn't stretched.
% - Set and use an API key which enables a much higher usage volume per day.
% Version 1.1 - 25/08/2011
persistent apiKey useTemp
if isnumeric(apiKey)
% first run, check if API key file exists
if exist('api_key.mat','file')
load api_key
else
apiKey = '';
end
end
if isempty(useTemp)
% first run, check we we have wrtie access to the temp folder
try
tempfilename = tempname;
fid = fopen(tempfilename, 'w');
if fid > 0
fclose(fid);
useTemp = true;
delete(tempfilename);
else
% Don't have write access to temp folder or it doesn't exist, fallback to current dir
useTemp = false;
end
catch
% in case tempname fails for some reason
useTemp = false;
end
end
hold on
% Default parametrs
axHandle = gca;
height = 640;
width = 640;
scale = 2;
maptype = 'roadmap';
alphaData = 1;
autoRefresh = 1;
figureResizeUpdate = 1;
autoAxis = 1;
showLabels = 1;
language = '';
markeridx = 1;
markerlist = {};
% Handle input arguments
if nargin >= 2
for idx = 1:2:length(varargin)
switch lower(varargin{idx})
case 'axis'
axHandle = varargin{idx+1};
case 'height'
height = varargin{idx+1};
case 'width'
width = varargin{idx+1};
case 'maptype'
maptype = varargin{idx+1};
case 'alpha'
alphaData = varargin{idx+1};
case 'refresh'
autoRefresh = varargin{idx+1};
case 'showlabels'
showLabels = varargin{idx+1};
case 'figureresizeupdate'
figureResizeUpdate = varargin{idx+1};
case 'language'
language = varargin{idx+1};
case 'marker'
markerlist{markeridx} = varargin{idx+1};
markeridx = markeridx + 1;
case 'autoaxis'
autoAxis = varargin{idx+1};
case 'apikey'
apiKey = varargin{idx+1}; % set new key
% save key to file
funcFile = which('plot_google_map.m');
pth = fileparts(funcFile);
keyFile = fullfile(pth,'api_key.mat');
save(keyFile,'apiKey')
otherwise
error(['Unrecognized variable: ' varargin{idx}])
end
end
end
if height > 640
height = 640;
end
if width > 640
width = 640;
end
% Store paramters in axis handle (for auto refresh callbacks)
ud = get(axHandle, 'UserData');
if isempty(ud)
% explicitly set as struct to avoid warnings
ud = struct;
end
ud.gmap_params = varargin;
set(axHandle, 'UserData', ud);
curAxis = axis(axHandle);
% Enforce Latitude constraints of EPSG:900913
if curAxis(3) < -85
curAxis(3) = -85;
end
if curAxis(4) > 85
curAxis(4) = 85;
end
% Enforce longitude constrains
if curAxis(1) < -180
curAxis(1) = -180;
end
if curAxis(1) > 180
curAxis(1) = 0;
end
if curAxis(2) > 180
curAxis(2) = 180;
end
if curAxis(2) < -180
curAxis(2) = 0;
end
if isequal(curAxis,[0 1 0 1]) % probably an empty figure
% display world map
curAxis = [-200 200 -85 85];
axis(curAxis)
end
if autoAxis
% adjust current axis limit to avoid strectched maps
[xExtent,yExtent] = latLonToMeters(curAxis(3:4), curAxis(1:2) );
xExtent = diff(xExtent); % just the size of the span
yExtent = diff(yExtent);
% get axes aspect ratio
drawnow
org_units = get(axHandle,'Units');
set(axHandle,'Units','Pixels')
ax_position = get(axHandle,'position');
set(axHandle,'Units',org_units)
aspect_ratio = ax_position(4) / ax_position(3);
if xExtent*aspect_ratio > yExtent
centerX = mean(curAxis(1:2));
centerY = mean(curAxis(3:4));
spanX = (curAxis(2)-curAxis(1))/2;
spanY = (curAxis(4)-curAxis(3))/2;
% enlarge the Y extent
spanY = spanY*xExtent*aspect_ratio/yExtent; % new span
if spanY > 85
spanX = spanX * 85 / spanY;
spanY = spanY * 85 / spanY;
end
curAxis(1) = centerX-spanX;
curAxis(2) = centerX+spanX;
curAxis(3) = centerY-spanY;
curAxis(4) = centerY+spanY;
elseif yExtent > xExtent*aspect_ratio
centerX = mean(curAxis(1:2));
centerY = mean(curAxis(3:4));
spanX = (curAxis(2)-curAxis(1))/2;
spanY = (curAxis(4)-curAxis(3))/2;
% enlarge the X extent
spanX = spanX*yExtent/(xExtent*aspect_ratio); % new span
if spanX > 180
spanY = spanY * 180 / spanX;
spanX = spanX * 180 / spanX;
end
curAxis(1) = centerX-spanX;
curAxis(2) = centerX+spanX;
curAxis(3) = centerY-spanY;
curAxis(4) = centerY+spanY;
end
% Enforce Latitude constraints of EPSG:900913
if curAxis(3) < -85
curAxis(3:4) = curAxis(3:4) + (-85 - curAxis(3));
end
if curAxis(4) > 85
curAxis(3:4) = curAxis(3:4) + (85 - curAxis(4));
end
axis(axHandle, curAxis); % update axis as quickly as possible, before downloading new image
drawnow
end
% Delete previous map from plot (if exists)
if nargout <= 1 % only if in plotting mode
curChildren = get(axHandle,'children');
map_objs = findobj(curChildren,'tag','gmap');
bd_callback = [];
for idx = 1:length(map_objs)
if ~isempty(get(map_objs(idx),'ButtonDownFcn'))
% copy callback properties from current map
bd_callback = get(map_objs(idx),'ButtonDownFcn');
end
end
delete(map_objs)
end
% Calculate zoom level for current axis limits
[xExtent,yExtent] = latLonToMeters(curAxis(3:4), curAxis(1:2) );
minResX = diff(xExtent) / width;
minResY = diff(yExtent) / height;
minRes = max([minResX minResY]);
tileSize = 256;
initialResolution = 2 * pi * 6378137 / tileSize; % 156543.03392804062 for tileSize 256 pixels
zoomlevel = floor(log2(initialResolution/minRes));
% Enforce valid zoom levels
if zoomlevel < 0
zoomlevel = 0;
end
if zoomlevel > 19
zoomlevel = 19;
end
% Calculate center coordinate in WGS1984
lat = (curAxis(3)+curAxis(4))/2;
lon = (curAxis(1)+curAxis(2))/2;
% Construct query URL
preamble = 'http://maps.googleapis.com/maps/api/staticmap';
location = ['?center=' num2str(lat,10) ',' num2str(lon,10)];
zoomStr = ['&zoom=' num2str(zoomlevel)];
sizeStr = ['&scale=' num2str(scale) '&size=' num2str(width) 'x' num2str(height)];
maptypeStr = ['&maptype=' maptype ];
if ~isempty(apiKey)
keyStr = ['&key=' apiKey];
else
keyStr = '';
end
markers = '&markers=';
for idx = 1:length(markerlist)
if idx < length(markerlist)
markers = [markers markerlist{idx} '%7C'];
else
markers = [markers markerlist{idx}];
end
end
if showLabels == 0
labelsStr = '&style=feature:all|element:labels|visibility:off';
else
labelsStr = '';
end
if ~isempty(language)
languageStr = ['&language=' language];
else
languageStr = '';
end
if ismember(maptype,{'satellite','hybrid'})
filename = 'tmp.jpg';
format = '&format=jpg';
convertNeeded = 0;
else
filename = 'tmp.png';
format = '&format=png';
convertNeeded = 1;
end
sensor = '&sensor=false';
url = [preamble location zoomStr sizeStr maptypeStr format markers labelsStr languageStr sensor keyStr];
% Get the image
if useTemp
filepath = fullfile(tempdir, filename);
else
filepath = filename;
end
try
urlwrite(url,filepath);
catch % error downloading map
warning(['Unable to download map form Google Servers.\n' ...
'Matlab error was: %s\n\n' ...
'Possible reasons: missing write permissions, no network connection, quota exceeded, or some other error.\n' ...
'Consider using an API key if quota problems persist.\n\n' ...
'To debug, try pasting the following URL in your browser, which may result in a more informative error:\n%s'], lasterr, url);
varargout{1} = [];
varargout{2} = [];
varargout{3} = [];
return
end
[M Mcolor] = imread(filepath);
M = cast(M,'double');
delete(filepath); % delete temp file
width = size(M,2);
height = size(M,1);
% Calculate a meshgrid of pixel coordinates in EPSG:900913
centerPixelY = round(height/2);
centerPixelX = round(width/2);
[centerX,centerY] = latLonToMeters(lat, lon ); % center coordinates in EPSG:900913
curResolution = initialResolution / 2^zoomlevel/scale; % meters/pixel (EPSG:900913)
xVec = centerX + ((1:width)-centerPixelX) * curResolution; % x vector
yVec = centerY + ((height:-1:1)-centerPixelY) * curResolution; % y vector
[xMesh,yMesh] = meshgrid(xVec,yVec); % construct meshgrid
% convert meshgrid to WGS1984
[lonMesh,latMesh] = metersToLatLon(xMesh,yMesh);
% We now want to convert the image from a colormap image with an uneven
% mesh grid, into an RGB truecolor image with a uniform grid.
% This would enable displaying it with IMAGE, instead of PCOLOR.
% Advantages are:
% 1) faster rendering
% 2) makes it possible to display together with other colormap annotations (PCOLOR, SCATTER etc.)
% Convert image from colormap type to RGB truecolor (if PNG is used)
if convertNeeded
imag = zeros(height,width,3);
for idx = 1:3
imag(:,:,idx) = reshape(Mcolor(M(:)+1+(idx-1)*size(Mcolor,1)),height,width);
end
else
imag = M/255;
end
% Next, project the data into a uniform WGS1984 grid
sizeFactor = 1; % factoring of new image
uniHeight = round(height*sizeFactor);
uniWidth = round(width*sizeFactor);
latVect = linspace(latMesh(1,1),latMesh(end,1),uniHeight);
lonVect = linspace(lonMesh(1,1),lonMesh(1,end),uniWidth);
[uniLonMesh,uniLatMesh] = meshgrid(lonVect,latVect);
uniImag = zeros(uniHeight,uniWidth,3);
% old version (projection using INTERP2)
% for idx = 1:3
% % 'nearest' method is the fastest. difference from other methods is neglible
% uniImag(:,:,idx) = interp2(lonMesh,latMesh,imag(:,:,idx),uniLonMesh,uniLatMesh,'nearest');
% end
uniImag = myTurboInterp2(lonMesh,latMesh,imag,uniLonMesh,uniLatMesh);
if nargout <= 1 % plot map
% display image
hold(axHandle, 'on');
cax = caxis;
h = image(lonVect,latVect,uniImag, 'Parent', axHandle);
caxis(cax); % Preserve caxis that is sometimes changed by the call to image()
set(axHandle,'YDir','Normal')
set(h,'tag','gmap')
set(h,'AlphaData',alphaData)
% add a dummy image to allow pan/zoom out to x2 of the image extent
h_tmp = image(lonVect([1 end]),latVect([1 end]),zeros(2),'Visible','off', 'Parent', axHandle);
set(h_tmp,'tag','gmap')
% older version (display without conversion to uniform grid)
% h =pcolor(lonMesh,latMesh,(M));
% colormap(Mcolor)
% caxis([0 255])
% warning off % to avoid strange rendering warnings
% shading flat
uistack(h,'bottom') % move map to bottom (so it doesn't hide previously drawn annotations)
axis(axHandle, curAxis) % restore original zoom
if nargout == 1
varargout{1} = h;
end
% if auto-refresh mode - override zoom callback to allow autumatic
% refresh of map upon zoom actions.
figHandle = axHandle;
while ~strcmpi(get(figHandle, 'Type'), 'figure')
% Recursively search for parent figure in case axes are in a panel
figHandle = get(figHandle, 'Parent');
end
zoomHandle = zoom(axHandle);
panHandle = pan(figHandle); % This isn't ideal, doesn't work for contained axis
if autoRefresh
set(zoomHandle,'ActionPostCallback',@update_google_map);
set(panHandle, 'ActionPostCallback', @update_google_map);
else % disable zoom override
set(zoomHandle,'ActionPostCallback',[]);
set(panHandle, 'ActionPostCallback',[]);
end
% set callback for figure resize function, to update extents if figure
% is streched.
if figureResizeUpdate &&isempty(get(figHandle, 'ResizeFcn'))
% set only if not already set by someone else
set(figHandle, 'ResizeFcn', @update_google_map_fig);
end
% set callback properties
set(h,'ButtonDownFcn',bd_callback);
else % don't plot, only return map
varargout{1} = lonVect;
varargout{2} = latVect;
varargout{3} = uniImag;
end
% Coordinate transformation functions
function [lon,lat] = metersToLatLon(x,y)
% Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum
originShift = 2 * pi * 6378137 / 2.0; % 20037508.342789244
lon = (x ./ originShift) * 180;
lat = (y ./ originShift) * 180;
lat = 180 / pi * (2 * atan( exp( lat * pi / 180)) - pi / 2);
function [x,y] = latLonToMeters(lat, lon )
% Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913"
originShift = 2 * pi * 6378137 / 2.0; % 20037508.342789244
x = lon * originShift / 180;
y = log(tan((90 + lat) * pi / 360 )) / (pi / 180);
y = y * originShift / 180;
function ZI = myTurboInterp2(X,Y,Z,XI,YI)
% An extremely fast nearest neighbour 2D interpolation, assuming both input
% and output grids consist only of squares, meaning:
% - uniform X for each column
% - uniform Y for each row
XI = XI(1,:);
X = X(1,:);
YI = YI(:,1);
Y = Y(:,1);
xiPos = nan*ones(size(XI));
xLen = length(X);
yiPos = nan*ones(size(YI));
yLen = length(Y);
% find x conversion
xPos = 1;
for idx = 1:length(xiPos)
if XI(idx) >= X(1) && XI(idx) <= X(end)
while xPos < xLen && X(xPos+1)<XI(idx)
xPos = xPos + 1;
end
diffs = abs(X(xPos:xPos+1)-XI(idx));
if diffs(1) < diffs(2)
xiPos(idx) = xPos;
else
xiPos(idx) = xPos + 1;
end
end
end
% find y conversion
yPos = 1;
for idx = 1:length(yiPos)
if YI(idx) <= Y(1) && YI(idx) >= Y(end)
while yPos < yLen && Y(yPos+1)>YI(idx)
yPos = yPos + 1;
end
diffs = abs(Y(yPos:yPos+1)-YI(idx));
if diffs(1) < diffs(2)
yiPos(idx) = yPos;
else
yiPos(idx) = yPos + 1;
end
end
end
ZI = Z(yiPos,xiPos,:);
function update_google_map(obj,evd)
% callback function for auto-refresh
drawnow;
try
axHandle = evd.Axes;
catch ex
% Event doesn't contain the correct axes. Panic!
axHandle = gca;
end
ud = get(axHandle, 'UserData');
if isfield(ud, 'gmap_params')
params = ud.gmap_params;
plot_google_map(params{:});
end
function update_google_map_fig(obj,evd)
% callback function for auto-refresh
drawnow;
axes_objs = findobj(get(gcf,'children'),'type','axes');
for idx = 1:length(axes_objs)
if ~isempty(findobj(get(axes_objs(idx),'children'),'tag','gmap'));
ud = get(axes_objs(idx), 'UserData');
if isfield(ud, 'gmap_params')
params = ud.gmap_params;
else
params = {};
end
% Add axes to inputs if needed
if ~sum(strcmpi(params, 'Axis'))
params = [params, {'Axis', axes_objs(idx)}];
end
plot_google_map(params{:});
end
end
Пример:
lat = [48.8708 51.5188 41.9260 40.4312 52.523 37.982];
lon = [2.4131 -0.1300 12.4951 -3.6788 13.415 23.715];
plot(lon,lat,'.r','MarkerSize',20)
plot_google_map
lon = [2.4131 -0.1300 12.4951 -3.6788 13.415 23.715];
plot(lon,lat,'.r','MarkerSize',20)
plot_google_map
Или
lat = [48.8708 51.5188 41.9260 40.4312 52.523 37.982];
lon = [2.4131 -0.1300 12.4951 -3.6788 13.415 23.715];
plot(lon,lat,'.r','MarkerSize',20)
plot_google_map('MapType', 'hybrid')
lon = [2.4131 -0.1300 12.4951 -3.6788 13.415 23.715];
plot(lon,lat,'.r','MarkerSize',20)
plot_google_map('MapType', 'hybrid')
Ссылки
- http://blogs.mathworks.com/pick/2012/05/04/plot-google-map/
- http://www.mathworks.com/matlabcentral/fileexchange/27627-zoharby-plot-google-map
- http://www.mathworks.com/matlabcentral/fileexchange/24113
- http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/
- http://developers.google.com/maps/documentation/staticmaps/
[ Хронологический вид ]Комментарии
Войдите, чтобы комментировать.