From 54f999150b128db9aa7f053582c3eacd62b0bee6 Mon Sep 17 00:00:00 2001 From: Marco Berzborn <marco.berzborn@akustik.rwth-aachen.de> Date: Tue, 6 Jun 2017 14:11:22 +0200 Subject: [PATCH] added map projection plot to itaCoordinates --- .../@itaCoordinates/plot_map_projection.m | 135 ++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 kernel/ClassStuff/@itaCoordinates/plot_map_projection.m diff --git a/kernel/ClassStuff/@itaCoordinates/plot_map_projection.m b/kernel/ClassStuff/@itaCoordinates/plot_map_projection.m new file mode 100644 index 00000000..44342506 --- /dev/null +++ b/kernel/ClassStuff/@itaCoordinates/plot_map_projection.m @@ -0,0 +1,135 @@ +function varargout = plot_map_projection(this, varargin) +%PLOT_MAP_PROJECTION - Plot a map projection of data on a sphere +% This function plots a map projection for a dataset on a full or +% partial spherical surface +% +% Syntax: +% plot_map_projection(this, data, options) +% +% Options (default): +% 'proj' ('wagner4') : Specifies the projection type (See mapping toolbox documentation for supported projections) +% 'phi0' (0) : Longitute angle phi0 = 0 -> x-axis at origin of the map +% 'shading' ('interp') : Specifies the shading algorithm +% 'db' (false) : db plot +% 'limits' ([]) : Colorbar limits +% +% Example: +% plot_map_projection(this, data, 'db', 'limits', [-50, 0]) +% +% See also: +% itaCoordinates, ita_sph_sampling, surf, scatter +% +% Reference page in Help browser +% <a href="matlab:doc plot_map_projection">doc plot_map_projection</a> + +% <ITA-Toolbox> +% This file is part of the ITA-Toolbox. Some rights reserved. +% You can find the license for this m-file in the license.txt file in the ITA-Toolbox folder. +% </ITA-Toolbox> + + +% Author: Marco Berzborn -- Email: marco.berzborn@akustik.rwth-aachen.de +% Created: 06-Jun-2017 + +sArgs = struct('pos1_data','double',... + 'phi0',0,... + 'proj','wagner4',... + 'db',false,... + 'shading','interp',... + 'limits',[],... + 'fgh',[],... + 'colormap',[]); +[data,sArgs] = ita_parse_arguments(sArgs,varargin); + +% convert to longitude, latitude coordinates required by the mapping toolbox +[azi,ele,~] = cart2sph(this.x,this.y,this.z); +lat = ele*180/pi; +% if phi0 = 0, x = 1 is the origin of the map at lat = 0, lon = 0 +lon = (azi-sArgs.phi0)*180/pi; +% longitude, latitude limits for the plot +latLim = [min(lat),max(lat)]; +lonLim = [min(lon),max(lon)]; + +if isempty(sArgs.fgh) + sArgs.fgh = figure(); +end +% project sampling coordinates +mstruct = defaultm(sArgs.proj); +mstruct.origin = [0,0]; +mstruct = defaultm(mstruct); +mstruct.maplonlimit = lonLim; +mstruct.maplatlimit = latLim; +mstruct.flinewidth = 1; +mstruct.flatlimit = latLim; +mstruct.flonlimit = lonLim; +[x,y] = mfwdtran(mstruct,lat,lon); + +LatTicks = 0:30:latLim(2); +LatTicks = [-fliplr(30:30:abs(latLim(1))),LatTicks(1:end)]; +LonTicks = 0:50:lonLim(2); +LonTicks = [-fliplr(50:50:abs(lonLim(1))),LonTicks(1:end)]; + +[latTicks,lonTicks] = meshgrid(latLim(1),LonTicks); +[XTicks,~] = mfwdtran(mstruct,latTicks,lonTicks); +[latTicks,lonTicks] = meshgrid(LatTicks,lonLim(2)); +[~,YTicks] = mfwdtran(mstruct,latTicks,lonTicks); + + +% calculate delaunay triangulation for the meshgrid +tri = delaunay(x,y); +if sArgs.db + data = 20*log10(abs(data)); +else + data = abs(data); +end +trisurf(tri,x,y,data); + +cb = colorbar; +cb.Label.String = 'Magnitude'; +if ~isempty(sArgs.limits) + caxis(sArgs.limits); +end + + +% use custom colormap if given +if ~isempty(sArgs.colormap) + colormap(sArgs.fgh,sArgs.colormap); +end +% set shading, will be set to interpolate by default +shading(sArgs.shading); + +% change view to top view +view([0,90]) + +% create axes corresponding to the chosen map projection +axh = axesm(sArgs.proj, ... + 'MapLatLim', latLim, ... + 'MapLonLim', lonLim, ... + 'Grid','on', ... + 'Frame','on', ... + 'LabelFormat','none', ... + 'PLineLocation',30, ... + 'MLineLocation',50, ... + 'Origin',[0,0]); +axh.YLim = latLim * pi/180; +axh.XLim = lonLim * pi/180; + +% set ticks +axh.YTick = YTicks; +axh.YTickLabels = LatTicks; +axh.XTick = XTicks; +axh.XTickLabels = LonTicks; + +axh.YTickLabel = strcat(axh.YTickLabel, '^\circ'); +axh.XTickLabel = strcat(axh.XTickLabel, '^\circ'); + +% set labels for x- and y-axes +axh.YLabel.String = 'Latitude'; +axh.XLabel.String = 'Longitude'; + +if nargout + varargout{1} = sArgs.fgh; +end + +%end function +end \ No newline at end of file -- GitLab