[Dart-dev] [4149] DART/trunk/models/wrf/matlab/map_wrf.m: better than nothing
nancy at ucar.edu
nancy at ucar.edu
Wed Nov 18 13:01:16 MST 2009
Revision: 4149
Author: thoar
Date: 2009-11-18 13:01:16 -0700 (Wed, 18 Nov 2009)
Log Message:
-----------
better than nothing
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 2009-11-17 21:46:19 UTC (rev 4148)
+++ DART/trunk/models/wrf/matlab/map_wrf.m 2009-11-18 20:01:16 UTC (rev 4149)
@@ -1,3 +1,24 @@
+function field = map_wrf(fname, varname, copystring, levelindx, timeindx )
+%% map_wrf simply creates a map of a WRF field.
+%
+% fname = 'Prior_Diag.nc';
+% varname = 'U_d01';
+% copystring = 'ensemble mean';
+% levelindx = 10;
+% timeindx = 1;
+%
+% Example:
+% map_wrf(fname, varname, copystring, levelindx, timeindx )
+
+if ( nargin == 3 )
+ levelindx = 1;
+ timeindx = 1;
+elseif (nargin == 4 )
+ timeindx = 1;
+elseif (nargin ~= 5 )
+ error('Wrong number of arguments ... must have 3,4, or 5')
+end
+
% Data Assimilation Research Testbed -- DART
% Copyright 2004-2007, Data Assimilation Research Section
% University Corporation for Atmospheric Research
@@ -9,167 +30,127 @@
% $Revision$
% $Date$
-% Select field to plot (U, V, W, GZ, T, MU, QV, QC, QR)
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist',fname)
+end
+if ( ~ nc_isvar(fname,varname) )
+ error('%s is not a variable in %s',varname,fname)
+end
-field_num = input('Input field type, 1=U, 2=V, 3=W, 4=GZ, 5=T, 6=MU, 7=QV, 8=QC, 9=QR: ');
+copyindex = get_copy_index(fname,copystring);
-map_proj = {'lambert', 'ups', 'mercator'};
+%% 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);
+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);
-% Get file name of true state file
-fname = 'True_State';
-xlon = getnc(fname, 'XLON');
-we = size(xlon, 2);
-xlat = getnc(fname, 'XLAT');
-sn = size(xlat, 1);
-level = getnc(fname, 'level');
-bt = size(level, 1);
-
-%--Read data
-nc = netcdf( 'True_State.nc' , 'read' ) ;
-stdlat1 = nc.TRUELAT1(:);
-stdlat2 = nc.TRUELAT2(:);
-cen_lat = nc.CEN_LAT(:);
-cen_lon = nc.CEN_LON(:);
-mp = nc.MAP_PROJ(:);
-
-close(nc)
-
minlat = min(xlat(:)); maxlat = max(xlat(:));
minlon = min(xlon(:)); maxlon = max(xlon(:));
-true_times = getnc(fname, 'time');
-num_true_times = size(true_times, 1)
+%% Set up map projection
+% may have to translate their projections to matlabs keywords
-stime = input('Initial time : ');
-ftime = input('End time : ');
+fprintf('map projection is %s\n',map_proj{mp})
-% Get level for free atmosphere fields
-if field_num == 6
- field_level = 1;
-else
- field_level = input('Input level: ');
-end
+%% Get the variable units, description, etc.
+% Determine variable to plot ...
-var_units = 'U (m/s)';
-var_name = 'U';
-maxlev = bt;
-iso = [0.5:1:5];
-if field_num > 1
-var_units = 'V (m/s)';
- var_name = 'V';
-end
-if field_num > 2
-var_units = 'W (m/s)';
- var_name = 'W';
- maxlev = bt + 1;
- iso = [0.01:0.01:0.1];
-end
-if field_num > 3
-var_units = 'GZ (m^2/s^2)';
- var_name = 'PH';
- iso = [50:50:300];
-end
-if field_num > 4
-var_units = 'T (K)';
- var_name = 'T';
- maxlev = bt;
- iso = [0.5:0.5:5];
-end
-if field_num > 5
-var_units = 'MU (Pa)';
- var_name = 'MU';
- maxlev = 1;
- iso = [100:100:600];
-end
-if field_num > 6
-var_units = 'QV (kg/kg)';
- var_name = 'QVAPOR';
- maxlev = bt;
- iso = [0.0001:0.0001:0.001];
-end
-if field_num > 7
-var_units = 'QC (kg/kg)';
- var_name = 'QCLOUD';
- iso = [0.00001:0.00001:0.0001];
-end
-if field_num > 8
-var_units = 'QR (kg/kg)';
- var_name = 'QRAIN';
- iso = [0.00001:0.00001:0.0001];
-end
+varinfo = nc_getvarinfo(fname,varname);
-scrsz = get(0,'ScreenSize');
-figure('Position',[1 scrsz(4)/2 0.8*scrsz(4) 0.8*scrsz(4)])
-
-m = ceil(sqrt(ftime-stime+1));
-
-pane = 1;
-
-for itime = stime:ftime
-
-plot_title = ['True state ' var_units ' Level: ' num2str(field_level) ' Time: ' num2str(itime)];
-
-if maxlev > 1
- stride = [1 1 1 1 1];
- corner_t = [itime -1 field_level -1 -1];
- end_point_t = [itime -1 field_level -1 -1];
-else
- stride = [1 1 1 1];
- corner_t = [itime -1 -1 -1];
- end_point_t = [itime -1 -1 -1];
+for i = 1:length(varinfo.Attribute)
+ attname = varinfo.Attribute(i).Name;
+ attvalu = varinfo.Attribute(i).Value;
+ switch lower(attname)
+ case 'name'
+ varinfo.name = attvalu;
+ case 'long_name'
+ varinfo.long_name = attvalu;
+ case 'units'
+ varinfo.units = attvalu;
+ end
end
-% Extract field
+% 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];
-stag_field = getnc(fname, var_name,corner_t,end_point_t,stride);
+for itime = timeindx
-field = zeros(sn,we);
-if field_num == 1
-for iy = 1:sn
-for ix = 1:we
- field(iy,ix) = (stag_field(iy,ix) + stag_field(iy,ix+1))/2.0;
-end
-end
-elseif field_num == 2
-for iy = 1:sn
-for ix = 1:we
- field(iy,ix) = (stag_field(iy,ix) + stag_field(iy+1,ix))/2.0;
-end
-end
-else
-field = stag_field;
-end
+ plot_title = sprintf('%s %s levelindex %d time %f',fname,varname,levelindx,times(timeindx));
-% Plot field
+ %% Determine the hyperslab indexing
+
+ myinfo.diagn_file = fname;
+ myinfo.copyindex = copyindex;
+ 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);
-subplot(m,m,pane);
+ 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
+
+ % Plot field
+
+ % subplot(m,m,pane);
+
+ axesm(map_proj{mp},'Origin',[0 cen_lon 0],'MapParallels',[stdlat1 stdlat2],...
+ 'MapLatLimit',[minlat maxlat],'MapLonLimit',[minlon maxlon]);
+ % framem;
+
+ states = shaperead('usastatelo', 'UseGeoCoords', true);
+ geoshow([states.Lat],[states.Lon],'Color','black');
+ % linem(states.Lon, states.Lat, 'color', [0 0 0]);
+
+ % axis( [-0.6 0.6 .2 1.2 ]) % This works pretty well for present CONUS domain
+
+ if min(field(:)) ~= max(field(:))
-axesm(map_proj{mp},'Origin',[0 cen_lon 0],'MapParallels',[stdlat1 stdlat2],...
- 'MapLatLimit',[minlat maxlat],'MapLonLimit',[minlon maxlon]); framem;
+ disp('contouring ...')
+
+ [C h]=contourm(xlat,xlon,field) ;
+ h = clabelm(C,h,'labelspacing',288); set(h,'Fontsize',9);
+
+ % Can get lat,lon lines via, e.g., contourm(xlat,xlon,xlat)
+ % [C h]=contourm(xlat,xlon,xlat,'k') ;
+ % h = clabelm(C,h,'labelspacing',288); set(h,'Fontsize',9);
+
+ title(plot_title,'Interpreter','none')
+ %colorbar('vert')
+ end
-plotm(coast,'color',[0 0 0]);
-plotm(usalo('statebvec'),'color',[0 0 0]);
-plotm(usalo('conusvec'),'color',[0 0 0]);
-
-% axis( [-0.6 0.6 .2 1.2 ]) % This works pretty well for present CONUS domain
-
-if min(field) ~= max(field)
-
-[C h]=contourm(xlat,xlon,field) ;
-h = clabelm(C,h,'labelspacing',288); set(h,'Fontsize',9);
-
-% Can get lat,lon lines via, e.g., contourm(xlat,xlon,xlat)
-[C h]=contourm(xlat,xlon,xlat,'k') ;
-h = clabelm(C,h,'labelspacing',288); set(h,'Fontsize',9);
-
-title(plot_title)
-%colorbar('vert')
-
end
-
-pane = pane + 1;
-
-end
-
-% Loop for another try
-%map_wrf;
More information about the Dart-dev
mailing list