Home > atmlab > geoplots > gmt > create_gmt_earth.m

create_gmt_earth

PURPOSE ^

CREATE_GMT_EARTH Main wrapper to gmt interface

SYNOPSIS ^

function commands = create_gmt_earth(in)

DESCRIPTION ^

 CREATE_GMT_EARTH Main wrapper to gmt interface

 See help gmt_plot for input details

 PURPOSE: Create cell strings that will be used in a system call.

 Created by Salomon Eliasson
 $Id: create_gmt_earth.m 8570 2013-08-10 18:36:48Z seliasson $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

create_gmt_earth.m

SOURCE CODE ^

0001 function commands = create_gmt_earth(in)
0002 % CREATE_GMT_EARTH Main wrapper to gmt interface
0003 %
0004 % See help gmt_plot for input details
0005 %
0006 % PURPOSE: Create cell strings that will be used in a system call.
0007 %
0008 % Created by Salomon Eliasson
0009 % $Id: create_gmt_earth.m 8570 2013-08-10 18:36:48Z seliasson $
0010 
0011 assert(exist('in','var')==1,'gmtlab:badInput',...
0012     'This is an internal function to be used by gmt_plot.m')
0013 
0014 in = setup_create_gmt_earth(in);
0015 
0016 F = in.filename;
0017 
0018 % Remove earlier .ps and gmt residules
0019 commands{1} = sprintf('rm -f %s.ps',F);
0020 %commands{end+1} = sprintf('rm -f .gmtcommands4 .gmtdefaults');
0021 
0022 % Set paper size
0023 commands{end+1} = 'gmtset PAPER_MEDIA A0+';
0024 
0025 % Annotations of axis
0026 commands{end+1} = sprintf('gmtset ANNOT_FONT_SIZE_PRIMARY %s',num2str(in.annot_font_size_primary));
0027 commands{end+1} = sprintf('gmtset BASEMAP_AXES %s',in.basemap_axis);
0028 commands{end+1} = sprintf('gmtset HEADER_OFFSET=%s',num2str(in.header_offset));
0029 commands{end+1} = sprintf('gmtset HEADER_FONT %s',num2str(in.header_font));
0030 commands{end+1} = sprintf('gmtset HEADER_FONT_SIZE %s',num2str(in.header_fontsize));
0031 commands{end+1} = sprintf('gmtset MEASURE_UNIT=%s',in.measure_unit);
0032 
0033 % other gmtset-commands you might have
0034 if isfield(in,'gmtset')
0035     for i = length(in.gmtset)
0036         commands{end+1} = in.gmtset{i}; %#ok<*AGROW>
0037     end
0038 end
0039 
0040 % PSBASEMAP (open PS-file)
0041 commands{end+1} = in.plot.basemap;
0042 
0043 % NEARNEIGHBOR
0044 if isfield(in,'nearneighbor') && ~isstruct(in.nearneighbor)
0045     commands{end+1} = in.nearneighbor;
0046     in.grdfile = [in.filename '.grd'];
0047 elseif ~in.gridded && isfield(in.plot,'grid')
0048     in.grdfile = [in.filename '.grd'];
0049     commands{end+1} = nearneighbor(in.grdfile,in.plot.grid);
0050 end
0051 
0052 % IMAGE (plot data)
0053 if isfield(in,'grdimage')
0054     commands{end+1} = in.grdimage;
0055 elseif ~in.nodata
0056     commands{end+1} = grdimage(F,in);
0057 end
0058 
0059 % LEGEND
0060 if isfield(in,'psscale')
0061     commands{end+1} = in.psscale;
0062 elseif ~isempty(in.plot.legend)
0063     commands{end+1} = psscale(F,in.plot.legend);
0064 end
0065 if isfield(in,'psscale_extra')
0066     commands{end+1} = in.psscale_extra;
0067 elseif ~isempty(in.plot.extralegend)
0068     commands{end+1} = in.plot.extralegend;
0069 end
0070 
0071 % COASTLINES
0072 if isfield(in,'pscoast')
0073     commands{end+1} = in.pscoast;
0074 elseif ~isempty(in.plot.coast)
0075     commands{end+1} = pscoast(F,in.plot.coast);
0076 end
0077 
0078 % LOCATION MARKERS
0079 if isfield(in,'locations')
0080     commands = pslocations(F,in,commands);
0081 end
0082 
0083 % CONTOURLINES
0084 if isfield(in,'grdcontour')
0085     commands{end+1} = in.grdcontour;
0086 elseif isfield(in.plot,'grdcontour')
0087     for i = 1:length(in.plot.grdcontour) %can have more than 1 set
0088         commands{end+1} = grdcontour(F,in.plot.grdcontour(i));
0089     end
0090 end
0091 
0092 % PSBOX & PSPOLY (draw boxes or shapes on map)
0093 if isfield(in,'psxy')
0094     commands{end+1} = in.psxy;
0095 elseif isfield(in,'psbox') || isfield(in.plot,'psbox')
0096     commands = psbox(F,in,commands);
0097 elseif isfield(in.plot,'pspoly')
0098     commands = pspoly(F,in,commands);
0099 end
0100 
0101 % PSTEXT
0102 if isfield(in,'pstext') && ~isstruct(in.pstext)
0103     commands{end+1} = in.pstext;
0104 elseif isfield(in.plot,'pstext')
0105     commands{end+1} = pstext(F,in.plot.pstext,in.plotPlacement);
0106 end
0107 
0108 % CLOSE PS + GRID
0109 commands{end+1} = in.plot.lastbasemap;
0110 
0111 %%%%%%%%%%%%%%%%%
0112 %  SUBFUNCTIONS
0113 %      |||||
0114 %      vvvvv
0115 function in = setup_create_gmt_earth(in)
0116 %% setup_create_gmt_earth %%
0117 
0118 in.plot.region = in.region;
0119 
0120 % PROJECTION
0121 if ~isfield(in,'proj')
0122     in.plot.proj = setup_projection(in);
0123 else in.plot.proj = in.proj;
0124 end
0125 
0126 in.plot.basemap = setupbasemap(in);
0127 
0128 % GRIDDING
0129 if ~in.gridded && ~in.nodata
0130     in.plot.grid = mkgrid(in);
0131 end
0132 
0133 % COLOR TABLE
0134 if ~in.nodata
0135     in.plot.color = create_colortables(in);
0136 end
0137 
0138 % COAST
0139 if ~isfield(in,'pscoast')
0140     in.plot.coast = mkcoast(in);
0141 end
0142 
0143 % PSBOX
0144 if isfield(in,'psbox')
0145     in.plot.psbox = setup_psboxpoly(in, 'psbox');
0146 end
0147 
0148 if isfield(in,'pspoly')
0149     in.plot.pspoly = setup_psboxpoly(in, 'pspoly');
0150 end
0151 
0152 % CONTOUR LINES
0153 if isfield(in,'contourline')
0154     in.plot.grdcontour = mkgrdcontour(in);
0155 end
0156 
0157 % PSTEXT
0158 if isfield(in,'pstext') && isstruct(in.pstext)
0159     in.plot.pstext = setup_pstext(in.pstext);
0160 end
0161 
0162 % LEGEND/S
0163 [in.plot.legend,in.plot.extralegend] = mklegend(in);
0164 
0165 % Close file with this
0166 in.plot.lastbasemap = last_psbasemap(in);
0167 
0168 function proj = setup_projection(in)
0169 %% setup_projection %%
0170 
0171 if isempty(in.center) || ~strcmp(in.plot.region(1:9),'-180/180/')
0172     logtext(1,'Data is not cyclic, choosing in.centre for projection from center longitude...\n')
0173     x = sscanf(in.plot.region,'%f/%f/%f/%f');
0174     C = (x(1)+x(2))/2;
0175 else C = in.center;
0176 end
0177 
0178 proj = sprintf('%s%g/%s',in.projection,C,num2str(in.map_width));
0179 
0180 function basemap = setupbasemap(in)
0181 %% setupbasemap
0182 
0183 if ~isempty(in.header)
0184     basemap = sprintf('psbasemap -R%s -J%s %s -G -P -K --COLOR_FOREGROUND=%s --COLOR_BACKGROUND=%s -B:."%s": > %s.ps',...
0185         in.plot.region,in.plot.proj,in.plotPlacement,in.color_foreground,in.color_background,gmt_unicode_converter(in.header),in.filename);
0186     
0187 else
0188     basemap = sprintf('psbasemap -R%s -J%s %s -G -P -K > %s.ps',...
0189         in.plot.region,in.plot.proj,in.plotPlacement,in.filename);
0190 end
0191 
0192 function a = mkgrid(in)
0193 %% mkgrid %%
0194 % If nothing is given the resolution and search radius is loosely based on the
0195 % density of the data points. nearneighbour uses the command: grdimage
0196 
0197 if isfield(in,'nearneighbor') && ~isstruct(in.nearneighbor)
0198     a='';
0199     return
0200 end
0201 a.ungriddedfile = in.ungriddedfile;
0202 a.plotPlacement = in.plotPlacement;
0203 
0204 % ----------------------
0205 % AREA
0206 
0207 % get number of 1deg boxes in region
0208 x            = sscanf(in.region,'%f/%f/%f/%f');
0209 oneDegRegion = ( x(2)-x(1) )*( x(4)-x(3) );
0210 
0211 
0212 % ---------------------
0213 % Maximum resolution
0214 
0215 % memGB = freeRAM()/1024; % How much RAM is available? TO LARGE
0216 % some reasonable maximum default
0217 % apparently it costs 1200bytes per point
0218 
0219 memGB = .1;
0220 maxpoints = memGB*1024^3  / 1200;
0221 maxres = oneDegRegion/maxpoints;
0222 
0223 % -----------------------------------------------------
0224 % Resolution inferred by data, or by user defined
0225 
0226 %data
0227 res = max([ abs(median(diff(in.lat)));abs(median(diff(in.lon)))]) ;
0228 
0229 %user
0230 if isfield(in,'nearneighbor') &&...
0231         isstruct(in.nearneighbor) &&...
0232         isfield(in.nearneighbor,'resolution')
0233     res = in.nearneighbor.resolution;
0234 end
0235 
0236 % check if res is acceptable
0237 if res < maxres
0238     logtext(atmlab('OUT'), ['Automatically picked the maximum resolution of: ',...
0239         '%.4fdeg instead of %.4fdeg defined by measurement density\n'],maxres,res)
0240     res = maxres;
0241 else
0242     logtext(atmlab('OUT'), 'Plot resolution: %.4f Deg\n',res)
0243 end
0244 a.increment = sprintf('%fm',60*res); % 60min x resolution
0245 
0246 % use increment to decide the search radius, unless given
0247 if ~isfield(in,'nearneighbor') || ~isfield(in.nearneighbor,'search')
0248     a.search = sprintf('%fm', 1.25*str2double(a.increment(1:end-1)));
0249 else
0250     a.search = in.nearneighbor.search;
0251 end
0252 
0253 function a = create_colortables(in)
0254 %% create color tables %%
0255 
0256 if isfield(in,'cptfile')
0257     a.cptfile = in.cptfile;
0258     return
0259 end
0260 
0261 % NOTE: remember that I am currently standing in a temporary directory
0262 
0263 % MAKE a custom colour table based on manually input color values
0264 if isfield(in,'colorrange')
0265     in.ctable = 'colorrange';
0266 end
0267 
0268 if isfield(in,'tickval')
0269     % make temporary tickval file
0270     fid=fopen('tickvalues.txt','w');
0271     if isfield(in,'tick_annotation_format')
0272         fprintf(fid,in.tick_annotation_format,in.tickval);
0273     else
0274         fprintf(fid,sprintf('%s\n',getAnnotFormat(in.tickval)),in.tickval);
0275     end
0276     fclose(fid);
0277 end
0278 
0279 % generate cpt files
0280 switch in.ctable
0281     case 'mypolar'
0282         if ~isfield(in,'tickval')
0283             tickval = ownctable(in);
0284         else tickval = in.tickval;
0285         end
0286         a.cptfile = makepolar(tickval,in);
0287     case 'colorrange'
0288         % make ctable to be input to makecpt -C%s
0289         if isfield(in.colorrange,'colors')
0290             cpt_from_colorrange(in.colorrange);
0291             a.cptfile = makecpt(in);
0292         elseif isfield(in.colorrange,'discrete_colors')
0293             a.cptfile = cpt_discrete_colors(in);
0294         end
0295         
0296     otherwise
0297         logtext(1,'Processing makecpt\n')
0298         if ~isfield(in,'makecpt')
0299             a.cptfile = makecpt(in);
0300         else
0301             exec_system_cmd(in.makecpt);
0302             tmp = splitstring(in.makecpt,'> ');
0303             a.cptfile = tmp{2};
0304         end
0305 end
0306 
0307 % extra legend box for missing values
0308 if isfield(in,'extra_legend')
0309     a.cptextra = extra_legend(in);
0310 end
0311 
0312 function a = mkcoast(in)
0313 %% mkcoast
0314 
0315 if ~isstruct(in.coast)
0316     a='';
0317     return
0318 end
0319 a.plotPlacement = in.plotPlacement;
0320 
0321 % Set resolution of coastal features
0322 rg=sscanf(in.plot.region,'%f/%f/%f/%f');
0323 ln=rg(2)-rg(1);
0324 lt=rg(4)-rg(3);
0325 if ~isfield(in.coast,'features')
0326     const = 65;%constant => features ~1000km for full map
0327     a.features = ln*lt/const;
0328 else a.features = in.coast.features;
0329 end
0330 if ~isfield(in.coast,'resolution')
0331     if ln < 20
0332         a.resolution = 'h';
0333     elseif ln < 45
0334         a.resolution = 'i';
0335     elseif ln < 90
0336         a.resolution = 'l';
0337     else a.resolution = 'c';
0338     end
0339 else a.resolution = in.coast.resolution;
0340 end
0341 
0342 if isfield(in.coast,'color')
0343     a.color = in.coast.color;
0344 end
0345 a.rivers = in.coast.rivers;
0346 a.width = in.coast.width;
0347 
0348 function a = mkgrdcontour(in)
0349 %% mkgrdcontour
0350 
0351 default.cptfile  = in.plot.color.cptfile;
0352 default.grdfile  = in.grdfile;
0353 
0354 for i = 1:length(in.contourline)
0355     a(i) = catstruct(default,in.contourline(i));
0356     a(i).plotPlacement = in.plotPlacement;
0357 end
0358 
0359 function out = setup_pstext(in)
0360 %% setup_pstext
0361 
0362 assert(all(isfield(in,{'text','lat','lon'})),['gmtlab:' mfilename ':BadInput'],...
0363     'in.pstext.{text,lat,lon} are the minimum input requirements for in.pstext')
0364 
0365 default.thick   = 20;      % text size in points
0366 default.angle   = 0;       % degrees counter-clockwise from horizontal
0367 default.fontnum = 1;       % sets the font type
0368 default.justify = 6;       % sets the alignment
0369 default.color   = '0/0/0'; % textcolor (black)
0370 for i = 1:length(in)
0371     out(i) = optargs_struct(in(i),default);
0372 end
0373 
0374 function [a,b] = mklegend(in)
0375 %% mklegend
0376 
0377 if isfield(in,'fieldname') && isstruct(in.legend) %if not field then nodata.
0378     a = in.legend;
0379     % get what is in the in-structure, i.e one level up from in.legend
0380     % structure where they really belong
0381     F = {'xunit','unit','cptfile','plotPlacement'};
0382     for f = F(ismember(F,fieldnames(in)))
0383         switch f{1}
0384             case {'unit','xunit'}
0385                 a.(f{1}) = gmt_unicode_converter(in.(f{1}));
0386             otherwise
0387                 a.(f{1}) = in.(f{1});
0388         end
0389     end
0390     % Also needed
0391     a.cptfile = in.plot.color.cptfile;
0392     if isfield(in,'extra_legend') && isstruct(in.extra_legend)
0393         if ~isfield(in.extra_legend,'position')
0394             [w,wu]=separate_integer_and_unit(in.legend.width); %thickness of legend bar
0395             [x,xu]=separate_integer_and_unit(in.legend.xpos);
0396             [y,yu]=separate_integer_and_unit(in.legend.ypos);
0397             [l,~]=separate_integer_and_unit(in.legend.length);
0398             if ~isequal(wu,xu,yu),warning(['atmlab:' mfilename],'miss-matching units will cause strange default placements of extra_legend');end
0399             if strcmp(in.legend.orientation,'v')
0400                 vec = [x, y-0.5*l-2*w, w, w];
0401                 lsize = sprintf('%g%s/%g%s/%g%s/%g%s',...
0402                     vec(1),xu,vec(2),wu,vec(3),wu,vec(3),wu);
0403             elseif strcmp(in.legend.orientation,'h')
0404                 vec = [x+0.5*l+2*w, y, w, w];
0405                 lsize = sprintf('%g%s/%g%s/%g%s/%g%s',...
0406                     vec(1),wu,vec(2),yu,vec(3),wu,vec(3),wu);
0407             end
0408         else
0409             lsize = in.extra_legend.position;
0410         end
0411         if isfield(in.extra_legend,'fontsize')
0412             fz = in.extra_legend.fontsize;
0413         else fz = in.annot_font_size_primary;
0414         end
0415         b = sprintf('psscale -D%s -Xa5 -Ya5 -Li -O -K -C%s --ANNOT_FONT_SIZE_PRIMARY=%s >> %s.ps',...
0416             lsize,in.plot.color.cptextra,fz,in.filename);
0417     end
0418     if isfield(in.legend,'tick_annotations')
0419         append_tickannotations(in.plot.color.cptfile,in.legend.tick_annotations);
0420     end
0421 end
0422 if ~exist('b','var'), b = struct([]);end
0423 if ~exist('a','var'), a = struct([]);end
0424 
0425 function tickval = ownctable(in)
0426 %% ownctable %%
0427 
0428 minn = in.datarange(1);maxx = in.datarange(2);
0429 if minn==maxx
0430     maxx=maxx+1;
0431 end
0432 tickval = minn:in.stepsize:maxx;
0433 
0434 function out = setup_psboxpoly(in, field)
0435 %% ps box
0436 
0437 el = size(in.(field),1);
0438 if isfield(in,[field 'color'])
0439     if length(in.([field 'color']))~=el
0440         error(['atmlab:' mfilename],'Need same number of colors as boxes')
0441     end
0442 else in.([field 'color'])(1:l) = {'k'}; %default is black
0443 end
0444 
0445 if isfield(in,([field 'thick']))
0446     if length(in.([field 'thick']))~=el
0447         error(['atmlab:' mfilename],'Need same number of sizes as boxes')
0448     end
0449 else out.thick(1:el)={10}; %default size=10
0450 end
0451 c = list_colors('colors',in.([field 'color']));
0452 colors = cell(length(c),1);
0453 for i=1:length(c)
0454     colors{i}=c{i}.*255; % instead of from 0:1
0455 end
0456 out.colors=colors;
0457 out.(field)=in.(field);
0458 
0459 function basemap = last_psbasemap(in)
0460 %% last_psbasemap
0461 
0462 % default.ticks = '60g60/30g30' (global)
0463 % in.plot.region has this format: 'lon1/lon2/lat1/lat2'
0464 rg=sscanf(in.plot.region,'%f/%f/%f/%f');
0465 lon=rg(2)-rg(1);
0466 lat=rg(4)-rg(3);
0467 
0468 % for nice step lengths
0469 
0470 % lons
0471 as   = [240 180 120 90 60 30 15];
0472 step = [60 40 30 20 15 10 5];
0473 a = step(~(as>lon));
0474 if isempty(a), a = lon/6; end
0475 
0476 %lats
0477 bs   = [120 60 90 40 30 20 15];
0478 step = [60 40 30 20 10 5 2.5];
0479 b = step(~(bs>lat));
0480 if isempty(b), b = lat/3; end
0481 
0482 default.ticks = sprintf('%gg%g/%gg%g',a(1),b(1),a(1)/2,b(1)/2);
0483 in = optargs_struct(in,default);
0484 
0485 basemap = sprintf('psbasemap -B%s',in.ticks);
0486 
0487 basemap = sprintf('%s -R -J %s -O >> %s.ps',basemap,in.plotPlacement,in.filename);

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005