[Dart-dev] [5839] DART/branches/development/observations: the matlab and scripts to generate synthetic obs_sequence
nancy at ucar.edu
nancy at ucar.edu
Wed Aug 8 09:40:18 MDT 2012
Revision: 5839
Author: nancy
Date: 2012-08-08 09:40:18 -0600 (Wed, 08 Aug 2012)
Log Message:
-----------
the matlab and scripts to generate synthetic obs_sequence
files evenly spaced on the globe.
Added Paths:
-----------
DART/branches/development/observations/even_sphere/
DART/branches/development/observations/even_sphere/README
DART/branches/development/observations/even_sphere/even_sphere.m
DART/branches/development/observations/even_sphere/even_sphere_dense.m
DART/branches/development/observations/even_sphere/even_sphere_pe2lyr.m
DART/branches/development/observations/even_sphere/run_fixed_network_daily.csh
DART/branches/development/observations/even_sphere/run_fixed_network_seq.csh
-------------- next part --------------
Added: DART/branches/development/observations/even_sphere/README
===================================================================
--- DART/branches/development/observations/even_sphere/README (rev 0)
+++ DART/branches/development/observations/even_sphere/README 2012-08-08 15:40:18 UTC (rev 5839)
@@ -0,0 +1,84 @@
+# DART software - Copyright 2004 - 2011 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
+#
+# <next few lines under version control, do not edit>
+# $URL$
+# $Id$
+# $Revision$
+# $Date$
+
+Generate a series of synthetic observations located at roughly
+evenly distributed locations on a sphere. At each location
+generate a vertical column of observations. This could mimic
+a radiosonde observing network, for example.
+
+
+
+This directory contains a MATLAB script that generates
+input for the 'create_obs_sequence' program. It takes
+a number of vertical levels and a total number of points,
+and generates a roughly evenly distributed set of observations
+across the entire globe. Note that the number of obs
+will be the number of points times the number of vertical
+levels.
+
+
+the process, end to end:
+
+MATLAB:
+
+edit even_sphere.m and set the number of levels, the
+number of profiles, the vertical coordinate type, etc.
+
+run it in MATLAB. it will make a plot (which you can
+save from the menu) and it will create a file 'even_create_input'.
+
+DART:
+
+build the following executables and have these files
+in the current directory:
+
+ ./create_obs_sequence
+ ./create_fixed_network_seq
+ input.nml
+
+(if these executables were compiled for a specific model,
+then if that model needs any other input files at startup
+time, they will need to be copied here as well.
+e.g. cam needs a caminput.nc and cam_phis.nc even though
+they will never be used.)
+
+1)
+run ./create_obs_sequence < even_create_input > /dev/null
+
+that makes a set_def.out file
+
+2)
+edit run_fixed_network_seq.csh to set the start/stop times
+
+run ./run_fixed_network_seq.csh which will call ./create_fixed_network_seq
+multiple times to make separate obs_seq files as output.
+this script is where you set the period between files.
+
+
+DETAILS on generating points evenly distributed on a sphere:
+
+this is the algorithm (i believe) that's being used:
+
+dlong := pi*(3-sqrt(5)) /* ~2.39996323 */
+dz := 2.0/N
+long := 0
+z := 1 - dz/2
+for k := 0 .. N-1
+ r := sqrt(1-z*z)
+ node[k] := (cos(long)*r, sin(long)*r, z)
+ z := z - dz
+ long := long + dlong
+
+(there's example code, in the 'My Golden Section Spiral' paragraph,
+in Python, here: http://www.xsi-blog.com/archives/115 )
+
+
+
+
Property changes on: DART/branches/development/observations/even_sphere/README
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/observations/even_sphere/even_sphere.m
===================================================================
--- DART/branches/development/observations/even_sphere/even_sphere.m (rev 0)
+++ DART/branches/development/observations/even_sphere/even_sphere.m 2012-08-08 15:40:18 UTC (rev 5839)
@@ -0,0 +1,161 @@
+% Generate approximately evenly distributed points on sphere using
+% Golden Section spiral algorithm
+% http://www.softimageblog.com/archives/115
+
+
+
+%% DART software - Copyright 2004 - 2011 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+close all; clear;
+
+% Generate obs_sequence input for this problem
+% Mandatory levels are : 1000 mb, 850 mb, 700 mb, 500 mb, 400 mb,
+% 300 mb, 200 mb, 150 mb, 100 mb, 50 mb, 30 mb, 20 mb, 10 mb, 7 mb, 5
+
+% For diagnostic test need a null data value and a null qc value
+diagnostic_test = false;
+
+% Input is in hectopascals
+levels = [1000 850 700 500 400 300 200 150 100 50 30 20 10 7 5];
+
+% Date information is overwritten by create_fixed_network_sequence
+year = 2008;
+month = 10;
+day = 1;
+hour = 0;
+
+% Number of roughly evenly distributed points in horizontal
+n = 600;
+x(1:n) = 0;
+y(1:n) = 0;
+z(1:n) = 0;
+
+% Total number of observations at single time is levels*n*3
+num_obs = size(levels, 2) * n * 3;
+
+inc = pi * (3 - sqrt(5));
+off = 2 / n;
+for k = 1:n
+ y(k) = (k-1) * off - 1 + (off / 2);
+ r = sqrt(1 - y(k)*y(k));
+ phi = (k-1) * inc;
+ x(k) = cos(phi) * r;
+ z(k) = sin(phi) * r;
+end
+
+% Now convert to latitude and longitude
+for k = 1:n
+ lon(k) = atan2(y(k), x(k)) + pi;
+ lat(k) = asin(z(k));
+ % Input is in degrees of lat and longitude
+ dlon(k) = rad2deg(lon(k));
+ dlat(k) = rad2deg(lat(k));
+end
+
+% Need to generate output for driving create_obs_sequence
+fid = fopen('even_create_input', 'w');
+
+% Output the total number of observations
+fprintf(fid, '%6i\n', num_obs);
+% 0 copies of data
+if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+else
+ fprintf(fid, '%2i\n', 0);
+end
+
+% 0 QC fields
+if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+else
+ fprintf(fid, '%2i\n', 0);
+end
+
+% Metadata for data copy
+if(diagnostic_test)
+ fprintf(fid, '"observations"\n');
+end
+% Metadata for qc
+if(diagnostic_test)
+ fprintf(fid, '"Quality Control"\n');
+end
+
+
+% Loop to create each observation
+for hloc = 1:n
+ for vloc = 1:size(levels, 2)
+ for field = 1:3
+ % 0 indicates that there is another observation;
+ fprintf(fid, '%2i\n', 0);
+ % Specify obs kind by string
+ if(field == 1)
+ fprintf(fid, 'RADIOSONDE_TEMPERATURE\n');
+ elseif(field == 2)
+ fprintf(fid, 'RADIOSONDE_U_WIND_COMPONENT\n');
+ elseif(field == 3)
+ fprintf(fid, 'RADIOSONDE_V_WIND_COMPONENT\n');
+ end
+ % Select pressure as the vertical coordinate
+ fprintf(fid, '%2i\n', 2);
+ % The vertical pressure level
+ fprintf(fid, '%5i\n', levels(vloc));
+ % Lon and lat in degrees next
+ fprintf(fid, '%6.2f\n', dlon(hloc));
+ fprintf(fid, '%6.2f\n', dlat(hloc));
+ % Now the date and time
+ fprintf(fid, '%5i %3i %3i %3i %2i %2i \n', year, month, day, hour, 0, 0);
+ % Finally, the error variance, 1 for temperature, 4 for winds
+ if(field == 1)
+ fprintf(fid, '%2i\n', 1);
+ else
+ fprintf(fid, '%2i\n', 4);
+ end
+
+ % Need value and qc for testing
+ if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+ fprintf(fid, '%2i\n', 0);
+ end
+ end
+ end
+end
+
+% File name
+fprintf(fid, 'set_def.out\n');
+
+% Done with output, close file
+fclose(fid);
+
+
+
+% Some plotting to visualize this observation set
+plot(lon, lat, '*');
+hold on
+
+% Plot an approximate CAM T85 grid for comparison
+% 256x128 points on A-grid
+% Plot latitude lines first
+del_lat = pi / 128;
+for i = 1:128
+ glat = del_lat * i - pi/2;
+ a = [0 2*pi];
+ b = [glat glat];
+ plot(a, b, 'k');
+end
+
+del_lon = 2*pi / 256;
+for i = 1:256
+ glon = del_lon*i;
+ a = [glon glon];
+ b = [-pi/2, pi/2];
+ plot(a, b, 'k');
+end
+
Property changes on: DART/branches/development/observations/even_sphere/even_sphere.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/observations/even_sphere/even_sphere_dense.m
===================================================================
--- DART/branches/development/observations/even_sphere/even_sphere_dense.m (rev 0)
+++ DART/branches/development/observations/even_sphere/even_sphere_dense.m 2012-08-08 15:40:18 UTC (rev 5839)
@@ -0,0 +1,168 @@
+% Generate approximately evenly distributed points on sphere using
+% Golden Section spiral algorithm
+% http://www.softimageblog.com/archives/115
+
+
+%% DART software - Copyright 2004 - 2011 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+close all; clear;
+
+% this differs from the original in that it quits at 150 mb (since
+% we rarely assimilate obs above this level), and it's more dense
+% than the original 600 points to generate more obs to assimilate.
+% i would like to do something about the 1000mb level since with
+% any topography it will be below the lowest pressure level in the
+% model. as a compromise, i am generating several extra levels near
+% the bottom which are NOT MANDATORY levels.
+
+% Generate obs_sequence input for this problem
+% levels used here are : 1000 mb, 950mb, 900mb, 850mb, 800 mb, 750 mb,
+% 700mb, 650mb, 600mb, 550mb, 500mb, 400 mb, 300 mb, 200 mb, 150 mb
+
+% For diagnostic test need a null data value and a null qc value
+diagnostic_test = false;
+
+% Input is in hectopascals
+levels = [1000 950 900 850 800 750 700 650 600 550 500 400 300 200 150];
+
+% Date information is overwritten by create_fixed_network_sequence
+year = 2008;
+month = 10;
+day = 1;
+hour = 0;
+
+% Number of roughly evenly distributed points in horizontal
+n = 2500;
+x(1:n) = 0;
+y(1:n) = 0;
+z(1:n) = 0;
+
+% Total number of observations at single time is levels*n*3
+num_obs = size(levels, 2) * n * 3;
+
+inc = pi * (3 - sqrt(5));
+off = 2 / n;
+for k = 1:n
+ y(k) = (k-1) * off - 1 + (off / 2);
+ r = sqrt(1 - y(k)*y(k));
+ phi = (k-1) * inc;
+ x(k) = cos(phi) * r;
+ z(k) = sin(phi) * r;
+end
+
+% Now convert to latitude and longitude
+for k = 1:n
+ lon(k) = atan2(y(k), x(k)) + pi;
+ lat(k) = asin(z(k));
+ % Input is in degrees of lat and longitude
+ dlon(k) = rad2deg(lon(k));
+ dlat(k) = rad2deg(lat(k));
+end
+
+% Need to generate output for driving create_obs_sequence
+fid = fopen('even_create_input', 'w');
+
+% Output the total number of observations
+fprintf(fid, '%6i\n', num_obs);
+% 0 copies of data
+if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+else
+ fprintf(fid, '%2i\n', 0);
+end
+
+% 0 QC fields
+if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+else
+ fprintf(fid, '%2i\n', 0);
+end
+
+% Metadata for data copy
+if(diagnostic_test)
+ fprintf(fid, '"observations"\n');
+end
+% Metadata for qc
+if(diagnostic_test)
+ fprintf(fid, '"Quality Control"\n');
+end
+
+
+% Loop to create each observation
+for hloc = 1:n
+ for vloc = 1:size(levels, 2)
+ for field = 1:3
+ % 0 indicates that there is another observation;
+ fprintf(fid, '%2i\n', 0);
+ % Specify obs kind by string
+ if(field == 1)
+ fprintf(fid, 'RADIOSONDE_TEMPERATURE\n');
+ elseif(field == 2)
+ fprintf(fid, 'RADIOSONDE_U_WIND_COMPONENT\n');
+ elseif(field == 3)
+ fprintf(fid, 'RADIOSONDE_V_WIND_COMPONENT\n');
+ end
+ % Select pressure as the vertical coordinate
+ fprintf(fid, '%2i\n', 2);
+ % The vertical pressure level
+ fprintf(fid, '%5i\n', levels(vloc));
+ % Lon and lat in degrees next
+ fprintf(fid, '%6.2f\n', dlon(hloc));
+ fprintf(fid, '%6.2f\n', dlat(hloc));
+ % Now the date and time
+ fprintf(fid, '%5i %3i %3i %3i %2i %2i \n', year, month, day, hour, 0, 0);
+ % Finally, the error variance, 1 for temperature, 4 for winds
+ if(field == 1)
+ fprintf(fid, '%2i\n', 1);
+ else
+ fprintf(fid, '%2i\n', 4);
+ end
+
+ % Need value and qc for testing
+ if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+ fprintf(fid, '%2i\n', 0);
+ end
+ end
+ end
+end
+
+% File name
+fprintf(fid, 'set_def.out\n');
+
+% Done with output, close file
+fclose(fid);
+
+
+
+% Some plotting to visualize this observation set
+plot(lon, lat, '*');
+hold on
+
+% Plot an approximate CAM T85 grid for comparison
+% 256x128 points on A-grid
+% Plot latitude lines first
+del_lat = pi / 128;
+for i = 1:128
+ glat = del_lat * i - pi/2;
+ a = [0 2*pi];
+ b = [glat glat];
+ plot(a, b, 'k');
+end
+
+del_lon = 2*pi / 256;
+for i = 1:256
+ glon = del_lon*i;
+ a = [glon glon];
+ b = [-pi/2, pi/2];
+ plot(a, b, 'k');
+end
+
Property changes on: DART/branches/development/observations/even_sphere/even_sphere_dense.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/observations/even_sphere/even_sphere_pe2lyr.m
===================================================================
--- DART/branches/development/observations/even_sphere/even_sphere_pe2lyr.m (rev 0)
+++ DART/branches/development/observations/even_sphere/even_sphere_pe2lyr.m 2012-08-08 15:40:18 UTC (rev 5839)
@@ -0,0 +1,162 @@
+% Generate approximately evenly distributed points on sphere using
+% Golden Section spiral algorithm
+% http://www.softimageblog.com/archives/115
+
+% customized for the pe2lyr model, which has 2 levels in
+% the vertical (and is global in lat/lon).
+
+%% DART software - Copyright 2004 - 2011 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+
+close all; clear;
+
+% Generate obs_sequence input for this problem
+% level 1 and 2
+
+% For diagnostic test need a null data value and a null qc value
+diagnostic_test = false;
+
+% Input is in level number
+levels = [1 2]
+
+% Date information is overwritten by create_fixed_network_sequence
+year = 2008;
+month = 10;
+day = 1;
+hour = 0;
+
+% Number of roughly evenly distributed points in horizontal
+n = 1600;
+x(1:n) = 0;
+y(1:n) = 0;
+z(1:n) = 0;
+
+% Total number of observations at single time is levels*n*3
+num_obs = size(levels, 2) * n * 3;
+
+inc = pi * (3 - sqrt(5));
+off = 2 / n;
+for k = 1:n
+ y(k) = (k-1) * off - 1 + (off / 2);
+ r = sqrt(1 - y(k)*y(k));
+ phi = (k-1) * inc;
+ x(k) = cos(phi) * r;
+ z(k) = sin(phi) * r;
+end
+
+% Now convert to latitude and longitude
+for k = 1:n
+ lon(k) = atan2(y(k), x(k)) + pi;
+ lat(k) = asin(z(k));
+ % Input is in degrees of lat and longitude
+ dlon(k) = rad2deg(lon(k));
+ dlat(k) = rad2deg(lat(k));
+end
+
+% Need to generate output for driving create_obs_sequence
+fid = fopen('even_create_input', 'w');
+
+% Output the total number of observations
+fprintf(fid, '%6i\n', num_obs);
+% 0 copies of data
+if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+else
+ fprintf(fid, '%2i\n', 0);
+end
+
+% 0 QC fields
+if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+else
+ fprintf(fid, '%2i\n', 0);
+end
+
+% Metadata for data copy
+if(diagnostic_test)
+ fprintf(fid, '"observations"\n');
+end
+% Metadata for qc
+if(diagnostic_test)
+ fprintf(fid, '"Quality Control"\n');
+end
+
+
+% Loop to create each observation
+for hloc = 1:n
+ for vloc = 1:size(levels, 2)
+ for field = 1:3
+ % 0 indicates that there is another observation;
+ fprintf(fid, '%2i\n', 0);
+ % Specify obs kind by string
+ if(field == 1)
+ fprintf(fid, 'U_WIND_COMPONENT\n');
+ elseif(field == 2)
+ fprintf(fid, 'V_WIND_COMPONENT\n');
+ elseif(field == 3)
+ fprintf(fid, 'GEOPOTENTIAL_HEIGHT\n');
+ end
+ % Select level as the vertical coordinate
+ fprintf(fid, '%2i\n', 1);
+ % The level number
+ fprintf(fid, '%2i\n', levels(vloc));
+ % Lon and lat in degrees next
+ fprintf(fid, '%6.2f\n', dlon(hloc));
+ fprintf(fid, '%6.2f\n', dlat(hloc));
+ % Now the date and time
+ fprintf(fid, '%5i %3i %3i %3i %2i %2i \n', year, month, day, hour, 0, 0);
+ % Finally, the error variance, 1 for temperature, 4 for winds
+ if(field == 1)
+ fprintf(fid, '%2i\n', 1);
+ else
+ fprintf(fid, '%2i\n', 4);
+ end
+
+ % Need value and qc for testing
+ if(diagnostic_test)
+ fprintf(fid, '%2i\n', 1);
+ fprintf(fid, '%2i\n', 0);
+ end
+ end
+ end
+end
+
+% File name
+fprintf(fid, 'set_def.out\n');
+
+% Done with output, close file
+fclose(fid);
+
+
+
+% Some plotting to visualize this observation set
+plot(lon, lat, '*');
+hold on
+
+% Plot an approximate CAM T85 grid for comparison
+% 256x128 points on A-grid
+% Plot latitude lines first
+del_lat = pi / 128;
+for i = 1:128
+ glat = del_lat * i - pi/2;
+ a = [0 2*pi];
+ b = [glat glat];
+ plot(a, b, 'k');
+end
+
+del_lon = 2*pi / 256;
+for i = 1:256
+ glon = del_lon*i;
+ a = [glon glon];
+ b = [-pi/2, pi/2];
+ plot(a, b, 'k');
+end
+
Property changes on: DART/branches/development/observations/even_sphere/even_sphere_pe2lyr.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/observations/even_sphere/run_fixed_network_daily.csh
===================================================================
--- DART/branches/development/observations/even_sphere/run_fixed_network_daily.csh (rev 0)
+++ DART/branches/development/observations/even_sphere/run_fixed_network_daily.csh 2012-08-08 15:40:18 UTC (rev 5839)
@@ -0,0 +1,37 @@
+#!/bin/csh
+
+# loop, calling ./create_fixed_network_seq to create separate
+# files for each time period. edit the values below to change
+# the dates and intervals.
+
+#
+# DART software - Copyright 2004 - 2011 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
+#
+# $Id$
+
+
+@ month = 8
+@ day = 1
+@ ndays = 31
+
+while($day <= $ndays)
+ echo '"set_def.out"' > create_fixed_input
+ echo 1 >> create_fixed_input
+ echo 2 >> create_fixed_input
+ echo "2008 ${month} ${day} 12 0 0" >> create_fixed_input
+ echo '0 43200' >> create_fixed_input
+
+ if($day < 10) then
+ echo obs_seq2008080${day} >> create_fixed_input
+ else
+ echo obs_seq200808${day} >> create_fixed_input
+ endif
+
+ ./create_fixed_network_seq < create_fixed_input
+ @ day++
+end
+
+exit 0
+
Property changes on: DART/branches/development/observations/even_sphere/run_fixed_network_daily.csh
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Added: DART/branches/development/observations/even_sphere/run_fixed_network_seq.csh
===================================================================
--- DART/branches/development/observations/even_sphere/run_fixed_network_seq.csh (rev 0)
+++ DART/branches/development/observations/even_sphere/run_fixed_network_seq.csh 2012-08-08 15:40:18 UTC (rev 5839)
@@ -0,0 +1,46 @@
+#!/bin/csh
+
+# loop, calling ./create_fixed_network_seq to create separate
+# files for each time period. edit the values below to change
+# the dates and intervals.
+
+# this one makes 2 files/day, 12 hours apart, single time per file.
+
+#
+# DART software - Copyright 2004 - 2011 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
+#
+# $Id$
+
+
+
+@ year = 2008
+@ month = 8
+@ day = 1
+@ ndays = 3
+@ nhours = 12
+
+
+while($day <= $ndays)
+ @ hour = 0
+ while ($hour < 24)
+ echo '"set_def.out"' > create_fixed_input
+ echo 1 >> create_fixed_input
+ echo 1 >> create_fixed_input
+ echo $year $month $day $hour 0 0 >> create_fixed_input
+ echo '0 0' >> create_fixed_input
+
+ set dstring = `printf %04d%02d%02d%02d $year $month $day $hour`
+ echo obs_seq${dstring} >> create_fixed_input
+
+ ./create_fixed_network_seq < create_fixed_input
+
+ @ hour += $nhours
+ end
+
+ @ day++
+end
+
+exit 0
+
Property changes on: DART/branches/development/observations/even_sphere/run_fixed_network_seq.csh
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
More information about the Dart-dev
mailing list