[Dart-dev] [4844] DART/trunk/models/wrf/matlab/map_wrf.m: Simplified pcolor image of any field in either a wrfinput_d0x or Prior/ Posterior

nancy at ucar.edu nancy at ucar.edu
Tue Apr 5 13:06:00 MDT 2011


Revision: 4844
Author:   thoar
Date:     2011-04-05 13:06:00 -0600 (Tue, 05 Apr 2011)
Log Message:
-----------
Simplified pcolor image of any field in either a wrfinput_d0x or Prior/Posterior
DART netCDF file. has examples on how to overlay coastlines, observation locations -
faithfully preserves the shape of the WRF domain.

Modified Paths:
--------------
    DART/trunk/models/wrf/matlab/map_wrf.m

-------------- next part --------------
Modified: DART/trunk/models/wrf/matlab/map_wrf.m
===================================================================
--- DART/trunk/models/wrf/matlab/map_wrf.m	2011-04-05 14:49:29 UTC (rev 4843)
+++ DART/trunk/models/wrf/matlab/map_wrf.m	2011-04-05 19:06:00 UTC (rev 4844)
@@ -1,17 +1,42 @@
-function h = map_wrf(fname, varname, copystring, levelindx, timeindx )
-%% map_wrf creates a map of a WRF field - please check what is supposed to
-%  happen with the staggering.
+function h = map_wrf(fname, varname, levelindx, timeindx, copystring )
+%% map_wrf creates an image of a WRF field without using the mapping toolbox.
 %
-% (rather slow) Example: 
-% fname      = 'Prior_Diag.nc';
+% Example using a DART-style diagnostic file, i.e.:
+%
+% fname      = '/glade/proj3/image/hliu/200812new/cwb_icbc/12_01/Prior_Diag.nc';
 % varname    = 'U_d01';
 % copystring = 'ensemble mean';
 % levelindx  = 10;
 % timeindx   = 1;
+% map_wrf(fname, varname, levelindx, timeindx, copystring);
+% worldmap;
+% axis off;
 %
-% map_wrf(fname, varname, copystring, levelindx, timeindx );
+%
+% Example using a WRF netCDF file:
+%
+% fname      = '/glade/proj3/image/hliu/ICBC_from_cwb/wrfinput_d01';
+% levelindx  = 10;
+% timeindx   = 1;
+% map_wrf(fname, 'U', levelindx, timeindx);
+% worldmap;   % superimpose some low-res coastlines 
+% axis off;
+%
+% layer on the locations of some observations:
+%
+% obsfile = 'obs_epoch_001.nc';
+% ObsTypeString = 'RADIOSONDE_U_WIND_COMPONENT';
+% region        = [0 360 -90 90 -Inf Inf];
+% CopyString    = 'observation';
+% QCString      = 'DART quality control';
+% verbose       = 1;   % anything > 0 == 'true'
+% obs = read_obs_netcdf(obsfile, ObsTypeString, region, CopyString, QCString, verbose);
+% hold on;
+% plot(obs.lons, obs.lats, 'kd', 'MarkerFaceColor','k')
+%
 
-%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+
+%% DART software - Copyright � 2004 - 2010 UCAR. This open source software is
 % provided by UCAR, "as is", without charge, subject to all terms of use at
 % http://www.image.ucar.edu/DAReS/DART/DART_download
 %
@@ -30,51 +55,21 @@
    error('Wrong number of arguments ... must have 3,4, or 5')
 end
 
+% Check to ensure the file and the desired variable exist
+
 if (exist(fname,'file') ~= 2) 
    error('%s does not exist',fname)
 end
 
-varexist(fname,{varname,'time','XLON_d01','XLAT_d01','level_d01', ...
-       'TRUELAT1','TRUELAT2', 'CEN_LAT', 'CEN_LON', 'MAP_PROJ'})
+varexist(fname,{varname})
 
-copyindex = get_copy_index(fname,copystring);
+%% Get the dimension information.
 
-%% Read map metadata
-times    = nc_varget(fname,     'time'); Ntimes = size(times, 1);
-xlon     = nc_varget(fname, 'XLON_d01');  we    = size( xlon, 2);
-xlat     = nc_varget(fname, 'XLAT_d01');  sn    = size( xlat, 1);
-level    = nc_varget(fname,'level_d01');  bt    = size(level, 1);
+xy          = GetWRFlatlon( fname, varname);
+levels      = GetWRFlevels( fname, varname);
+times       = GetWRFtimes(  fname, varname);
+timestrings = datestr(times,31);
 
-calendar    = nc_attget(fname, 'time','calendar');
-timeunits   = nc_attget(fname, 'time','units');
-timebase    = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
-timeorigin  = datenum(timebase(1),timebase(2),timebase(3));
-timestrings = datestr(times + timeorigin, 0);
-
-stdlat1  = nc_varget(fname, 'TRUELAT1');
-stdlat2  = nc_varget(fname, 'TRUELAT2');
-cen_lat  = nc_varget(fname,  'CEN_LAT');
-cen_lon  = nc_varget(fname,  'CEN_LON');
-mp       = nc_varget(fname, 'MAP_PROJ')+1;  % convert [0,N-1] to [1,N]
-map_proj = nc_attget(fname, 'MAP_PROJ','units');
-inds     = find((map_proj == ',') | (map_proj == '='));
-map_proj(inds) = ' ';
-[a b] = strread(map_proj,'%s %s','delimiter',' ');
-map_proj = lower(b);
-
-%% Check to make sure longitudes are [0, 360]
-inds = xlon < 0;
-xlon(inds) = xlon(inds) + 360;
-minlat = min(xlat(:)); maxlat = max(xlat(:));
-minlon = min(xlon(:)); maxlon = max(xlon(:));
-hunk.xlat = xlat;
-hunk.xlon = xlon;
-
-%% Set up map projection
-%  may have to translate their projections to matlabs keywords
-
-fprintf('map projection is %s\n',map_proj{mp})
-
 %% Get the variable units, description, etc.
 % Determine variable to plot ...
 
@@ -93,15 +88,14 @@
    end
 end
 
-% iso = [0.5:1:5];
-% iso = [0.01:0.01:0.1];
-% iso = [50:50:300];
-% iso = [0.5:0.5:5];
-% iso = [100:100:600];
-% iso = [0.0001:0.0001:0.001];
-% iso = [0.00001:0.00001:0.0001];
-% iso = [0.00001:0.00001:0.0001];
+copydim = find(strncmp('copy',varinfo.Dimension,length('copy')));
 
+if (isempty(copydim))
+   copyindex = NaN;
+else
+   copyindex = get_copy_index(fname,copystring);
+end
+
 for itime = timeindx
 
    plot_title = {fname, ...
@@ -110,72 +104,23 @@
    %% Determine the hyperslab indexing
    
    myinfo.diagn_file = fname;
-   myinfo.copyindex  = copyindex;
+   if (isfinite(copyindex)), myinfo.copyindex  = copyindex; end
    myinfo.levelindex = levelindx;
    myinfo.timeindex  = timeindx;
    [start, count]    = GetNCindices(myinfo, 'diagn', varname);
    
    % Extract field
    
-   stag_field = nc_varget(fname, varname, start, count);
-   field      = zeros(sn,we);
-
-   switch lower(varname)
-     case {'u','u10_d01','u_d01'}
-        for iy = 1:sn
-        for ix = 1:we
-           field(iy,ix) = (stag_field(iy,ix) + stag_field(iy,ix+1))/2.0;
-        end
-        end
-     case {'v','v10_d01','v_d01'}
-        for iy = 1:sn
-        for ix = 1:we
-           field(iy,ix) = (stag_field(iy,ix) + stag_field(iy+1,ix))/2.0;
-        end
-        end
-     otherwise
-        field = stag_field;
-   end
+   datmat = double(nc_varget(fname, varname, start, count));
    
-   % Plot field
+   clf; 
+   h = pcolor(xy.lonmat, xy.latmat, datmat);
+   set(h,'LineStyle','none');
    
-   % subplot(m,m,pane);
+   title(plot_title,'Interpreter','none')
+   h2 = colorbar('vert');
+   set(get(h2,'YLabel'),'String',varinfo.units)
 
-   switch lower(map_proj{mp})
-   case 'mercator'
-      axesm(map_proj{mp}, ...
-            'MapParallels',[stdlat1], ...
-            'MapLatLimit',[minlat maxlat], ...
-            'MapLonLimit',[minlon maxlon]);
-      h = gridm('on');
-      % mlabel on;
-      % plabel on;
-   
-   otherwise
-      fprintf('projection %s is not yet supported\n',map_proj{mp})
-   end
-   
-%  states = shaperead('usastatelo', 'UseGeoCoords', true);
-%  geoshow([states.Lat],[states.Lon],'Color','black');
-
-   landareas = shaperead('landareas', 'UseGeoCoords', true);
-   geoshow([landareas.Lat],[landareas.Lon],'Color','black');
-
-   if min(field(:)) ~= max(field(:))
-
-%     disp('contouring (slow) ...')
-%     [C h]=contourm(xlat,xlon,field) ;
-%     h = clabelm(C,h,'labelspacing',288);  set(h,'Fontsize',9);
-   
-      h1 = surfm(xlat,xlon,field);
-   
-      title(plot_title,'Interpreter','none')
-      h2 = colorbar('vert');
-      set(get(h2,'YLabel'),'String',varinfo.units)
-   else
-      error('field is all uniform value')
-   end
-
 end
 
 
@@ -190,7 +135,7 @@
 for i = 1:nvars
    gotone(i) = nc_isvar(filename,varnames{i});
    if ( ~ gotone(i) )
-      fprintf('\n%s is not a variable in %s\n',varnames{i},filename)
+      warning('\n%s is not a variable in %s\n',varnames{i},filename)
    end
 end
 
@@ -198,3 +143,61 @@
    error('missing required variable ... exiting')
 end
 
+
+
+function xy = GetWRFlatlon(fname, varname);
+%% Each of the WRF variables has a 'coordinate' attribute signifying which
+% of the 6 possible lat/lon variables is appropriate.
+
+coordinates{1} = sscanf(nc_attget(fname,varname,'coordinates'),'%s %*s');
+coordinates{2} = sscanf(nc_attget(fname,varname,'coordinates'),'%*s %s');
+
+latcoord       = find(strncmp('XLAT', coordinates, length('XLAT')) > 0);
+loncoord       = find(strncmp('XLON', coordinates, length('XLON')) > 0);
+xy.latmat      = nc_varget(fname, coordinates{latcoord});
+xy.lonmat      = nc_varget(fname, coordinates{loncoord});
+xy.latunits    = nc_attget(fname, coordinates{latcoord},'units');
+xy.lonunits    = nc_attget(fname, coordinates{latcoord},'units');
+
+inds = (xy.lonmat < 0); % Convert to 0,360 to minimize dateline probs.
+xy.lonmat(inds) = xy.lonmat(inds) + 360.0;
+
+
+
+function levels = GetWRFlevels(fname, varname);
+%% Return the appropriate vertical indices.
+
+varinfo = nc_getvarinfo(fname,varname);
+lvlind  = find(strncmp('bottom_top',varinfo.Dimension,length('bottom_top')));
+nlevels = varinfo.Size(lvlind);
+levels  = 1:nlevels;
+
+
+
+function times = GetWRFtimes(fname, varname);
+%% Return the appropriate time coordinates.
+%  DART files use 'time' as a real array
+%  WRF  files use 'Time' as a character string
+
+varinfo = nc_getvarinfo(fname,varname);
+timeind = find(strncmp('time',varinfo.Dimension,length('time')));
+
+if (isempty(timeind)) % WRF flavor Time variable.
+    timeind = find(strncmp('Time',varinfo.Dimension,length('Time')));
+    if (isempty(timeind)), error('%s has no time information. Dying.'); end
+    timestr  = nc_varget(fname, 'Times');  % Note the plural ... seeth.
+    timebase = sscanf(timestr,'%d%*c%d%*c%d%*c%d%*c%d%*c%d'); % YYYY MM DD HH MM SS
+    times    = datenum(timebase');
+
+else
+    times       = nc_varget(fname, 'time');
+    calendar    = nc_attget(fname, 'time','calendar');
+    timeunits   = nc_attget(fname, 'time','units');
+    timebase    = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+    timeorigin  = datenum(timebase(1),timebase(2),timebase(3));
+    times       = times + timeorigin;
+ 
+end
+
+
+


More information about the Dart-dev mailing list