[cig-commits] r18610 - in seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot: gmt matlab matlab/matlab_scripts
carltape at geodynamics.org
carltape at geodynamics.org
Mon Jun 13 15:18:31 PDT 2011
Author: carltape
Date: 2011-06-13 15:18:30 -0700 (Mon, 13 Jun 2011)
New Revision: 18610
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m
Log:
updated plotting script for horizontal cross sections (ADJOINT_TOMO)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2011-06-13 21:59:19 UTC (rev 18609)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2011-06-13 22:18:30 UTC (rev 18610)
@@ -4,7 +4,7 @@
#
# plot_horz_models.pl
# Carl Tape
-# 07-Dec-2008
+# 07-Dec-2008 (last used 13-June-2011)
#
# This script reads in a horizontal cross-section data file with five columns:
# lon, lat, model-1, model-2, ln(model-2 / model-1)
@@ -43,7 +43,7 @@
$iletter = 0; # label for publication figures
$letter = "C";
$iextra = 0; # extra information on figures
-$irun = 1; # default 2; 1 for m01/m00
+$irun = 2; # default 2; 1 for m01/m00
$imask = 1;
$isingle = 0; # single figure showing only final model
@@ -58,7 +58,7 @@
$dirdat = "INPUT/horz_${stirun}";
# KEY: file showing the cuts and the plotting range for the values
-$fcuts = "${dirdat}/horz_${stirun}_xc_cuts_mod";
+$fcuts = "${dirdat}/horz_xc_cuts_mod";
if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
# load values and check that the files exist
@@ -70,7 +70,7 @@
$vsnorms[$p-1] = $vs_norm;
$vbnorms[$p-1] = $vb_norm;
$zcuts[$p-1] = -$zcut;
- $dfiles[$p-1] = "${dirdat}/horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}_mask${imask}.dat";
+ $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}_mask${imask}.dat";
#if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
}
print "\nNumber of cuts is $nump\n";
@@ -103,7 +103,7 @@
#print "$dtitle\n"; die("TESTING");
#---------------------------------------------------------
-# DATA FILES
+# DATA FILES (MODIFY FOR EACH USER)
$dir0 = "/home/carltape/gmt";
$dir = "socal_2005";
@@ -204,8 +204,8 @@
#$coast_info = "-A5000 -Dl -W1.0p -S220/255/255 -C220/255/255 -G200"; # A : smallest feature plotted, in km^2; D : resolution
$coast_info = "${coast_res} -W0.75p -Na/0.75p";
-$fault_info = "-M -W0.75p,255/0/0";
-$fault_infoK = "-M -W0.75p,0/0/0";
+$fault_info = "-m -W0.75p,255/0/0";
+$fault_infoK = "-m -W0.75p,0/0/0";
$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
@@ -226,7 +226,7 @@
$cshfile = "plot_horz_models.csh";
open(CSH,">$cshfile");
-print CSH "gmtset PAPER_MEDIA letter BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tick LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1 FRAME_PEN $fpen TICK_PEN $tpen \n";
+print CSH "gmtset PAPER_MEDIA letter MEASURE_UNIT inch BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tick LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1 FRAME_PEN $fpen TICK_PEN $tpen \n";
@shifts = ($origin,"-X$dX","-X$dX");
if($ilonlab==1) {@iBs = (2,14,2);} else {@iBs = (15,9,15);}
@@ -234,7 +234,7 @@
for ($p = $pmin; $p <= $pmax; $p ++ ) {
$stip = sprintf("%3.3i",$p);
- $fname0 = "horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}";
+ $fname0 = "horz_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}";
$fname = $fname0;
$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
@@ -251,7 +251,7 @@
# datafile -- COLUMNS: lon, lat, m00, m16
$dfile = $dfiles[$p-1];
if (not -f $dfile) {die("Check if dfile $dfile exist or not\n");}
- print "Data file is $dfile\n";
+ print "DATA FILE: $dfile\n";
if($ivs == 1) {$cnorm = $vsnorm;} else {$cnorm = $vbnorm;}
$cpert1 = $vpert;
@@ -289,8 +289,9 @@
#$gradfile0 = "/home/datalib/Topography/GLOBAL/ETOPO1/ETOPO1_Ice_g.grad";
#if (not -f $gradfile0) {die("Check if gradfile $gradfile0 exist or not\n");}
+ # LOOP over initial model, final model, and logarithmic difference
for ($k = 1; $k <= 3; $k ++ ) {
- #for ($k = 1; $k <= 1; $k ++ ) {
+ #for ($k = 1; $k <= 1; $k ++ ) {
$shift = $shifts[$k-1];
$iB = $iBs[$k-1];
@@ -302,13 +303,14 @@
} else {
print CSH "psbasemap $J $R $B -K -O -V $shift >> $psfile\n";
}
- print CSH "psbasemap $J $R $B $cgray -K -O -V >> $psfile\n";
+ print CSH "psbasemap $J $R $B $cgray -K -O -V >> $psfile\n"; # mask color
if ($icolor==1) {
if ($k == 1) {
##print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
#print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
#print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$3/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$3/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grep -v NaN $dfile | awk '{print \$1,\$2,\$3/1000}' | xyz2grd -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1b $J -Q -K -O -V >> $psfile\n";
#print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
@@ -323,7 +325,7 @@
##print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
#print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
#print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
- print CSH "awk '{print \$1,\$2,\$4/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grep -v NaN $dfile | awk '{print \$1,\$2,\$4/1000}' | xyz2grd -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1b $J -Q -K -O -V >> $psfile\n";
#print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
@@ -334,7 +336,7 @@
print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 $fsize2 0 $fontno LM $slab2\nEOF\n";
- } else {
+ } elsif ($k == 3) {
#print CSH "awk '{print \$1,\$2,log(\$4/\$3)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
#print CSH "awk '{print \$1,\$2,\$5}' $dfile | nearneighbor -G$grdfile $R $interp\n";
print CSH "awk '{print \$1,\$2,\$5}' $dfile | xyz2grd -G$grdfile $R $interp\n";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m 2011-06-13 21:59:19 UTC (rev 18609)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m 2011-06-13 22:18:30 UTC (rev 18610)
@@ -1,29 +1,32 @@
%
-% function
+% function ax_out = axes_utm2ll(ax_in,s_zone,i_type)
% Carl Tape, 01-Aug-2007
-% printed xxx
%
% This function converts between a bounding region in lat-lon and in UTM.
%
% EXAMPLE:
-
-% ax_utm = axes_utm2ll([-121.6 -114.7 32.2 36.8],11,0)
-%
+% ax_utm = axes_utm2ll([-121.6 -114.7 32.2 36.8],'11S',0)
%
% calls utm2ll.m
% called by xxx
%
-function ax_out = axes_utm2ll(ax_in,i_zone,i_type)
+function ax_out = axes_utm2ll(ax_in,s_zone,i_type,ellipsoid)
xmin0 = ax_in(1);
xmax0 = ax_in(2);
ymin0 = ax_in(3);
ymax0 = ax_in(4);
-[xmin,ymin] = utm2ll(xmin0,ymin0,i_zone,i_type);
-[xmax,ymax] = utm2ll(xmax0,ymax0,i_zone,i_type);
+% if no ellipsoid is given then use Matlab default for that zone
+if nargin == 3
+ disp('no ellipsoid is provided as input');
+ [ellipsoid,estr] = utmgeoid(s_zone)
+end
+[xmin,ymin] = utm2ll(xmin0,ymin0,s_zone,i_type,ellipsoid);
+[xmax,ymax] = utm2ll(xmax0,ymax0,s_zone,i_type,ellipsoid);
+
ax_out = [xmin xmax ymin ymax];
if 0==1
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m 2011-06-13 21:59:19 UTC (rev 18609)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m 2011-06-13 22:18:30 UTC (rev 18610)
@@ -1,171 +1,93 @@
%
-% function [xo,yo,i_zone] = utm2ll(xi,yi,i_zone,i_type)
-% Carl Tape, 22-April-2004
-% printed xxx
+% function [x,y,estr] = utm2ll(xi,yi,s_zone,i_type)
%
-% This function converts between lat-lon and UTM coordinates.
-% Note: input is NOT in vector form (one value at a time).
+% Black-box Matlab code to convert between utm and lon-lat coordinates.
+% NOTE: Requires Matlab Mapping Toolbox.
%
-% FUNCTION INPUT:
-% xi longitude or x-UTM
-% yi latitude or y-UTM
-% i_type lat-lon --> UTM (=0) or UTM --> lat-lon (=1)
-% i_zone UTM zone (if i_type=1)
+% To find a UTM zone for a specific point, use utmzone(lat,lon)
+% UTM zones: http://www.dmap.co.uk/utmworld.htm
%
+% Ellipsoids available: "help almanac"
+% everest 1830 Everest ellipsoid
+% bessel 1841 Bessel ellipsoid
+% airy 1849 Airy ellipsoid
+% clarke66 1866 Clarke ellipsoid
+% clarke80 1880 Clarke ellipsoid
+% international 1924 International ellipsoid
+% krasovsky 1940 Krasovsky ellipsoid
+% wgs60 1960 World Geodetic System ellipsoid
+% iau65 1965 International Astronomical Union ellipsoid
+% wgs66 1966 World Geodetic System ellipsoid
+% iau68 1968 International Astronomical Union ellipsoid
+% wgs72 1972 World Geodetic System ellipsoid
+% grs80 1980 Geodetic Reference System ellipsoid
+% wgs84 1984 World Geodetic System ellipsoid
%
-%****************************************************************
-%**
-%** FILE NAME: utmtoll
-%**
-%** DATE WRITTEN: 7/22/93
-%**
-%** PROGRAMMER: Scott Hensley
-%**
-%** FUNCTIONAL DESCRIPTION: This routine converts between lat
-%** lon and utm coordinates for a datum determined from the input
-%** a and e2.
-%**
-%** ROUTINES CALLED:none
-%**
-%** NOTES: none
-%**
-%** UPDATE LOG:
-%** MATLAB version, Yuri Fialko 11/25/99
-%****************************************************************
+% Mapping Toolbox ellipsoid representations are two-element vectors,
+% called ellipsoid vectors. The ellipsoid vector has the form
+% [semimajor_axis eccentricity]. The semimajor axis can be in any
+% unit of distance; the choice of units typically drives the units
+% used for distance outputs in the toolbox functions.
+% Meters, kilometers, or Earth radii (i.e., a unit sphere) are most
+% frequently used. See Functions that Define Ellipsoid Vectors for details.
%
-% INPUT VARIABLES:
-% real*8 xi,yi !input coordinates
-% integer i_zone !UTM zone (needed if i_type=2)
-% integer i_type ! [1 = lat,lon to utm], [2 = utm to lat,lon]
-% old stuff:
-% real*8 r_a !ellispoid semi-major axis
-% real*8 r_e2 !ellipsoid eccentricity squared
-% real*8 r_v(2) !geocentric vector (meters)
-% real*8 r_lat !latitude (deg -90 to 90)
-% real*8 r_lon !longitude (deg -180 to 180)
-% character*1 a_grid !UTM North-South grid
-%
-% OUTPUT VARIABLES:
-% real*8 xo,yo !output coordinates
+% calls xxx
+% called by test_utm.m
%
-% LOCAL VARIABLES:
-% integer i_ft,i_gi
-% real*8 pi,r_dtor
-% real*8 r_ep2,r_k0,r_k
-% real*8 r_fe,r_fn(2)
-% real*8 r_e4,r_e6,r_n,r_t,r_t2,r_c,r_c2,r_ba
-% real*8 r_a2,r_a3,r_a4,r_a5,r_a6
-% real*8 r_d,r_d2,r_d3,r_d4,r_d5,r_d6
-% real*8 r_lon0,r_lat1,r_m,r_m0,r_mu,r_lat0
-% real*8 r_et,r_e1,r_e12,r_e13,r_e14,r_r,r_vu(2)
-% character*1 a_griddes(20)
-%
-%****************************************************************
-function [xo,yo,i_zone] = utm2ll(xi,yi,i_zone,i_type)
+function [x,y,estr] = utm2ll(xi,yi,s_zone,i_type,ellipsoid)
-if and(i_type ~= 1, i_type ~= 0), error('i_type must be 1 or 0'); end
+utmstruct = defaultm('utm');
+utmstruct.zone = s_zone; % e.g., '11S'
-% DATA STATEMENTS
+% if no ellipsoid is given then use Matlab default ellipsoid for that zone
+if nargin == 4
+ disp('no ellipsoid is provided as input -- using first listed one in utmgeoid.m');
+ % NOTE: utmgeoid may return MULTIPLE geoids for a given zone
+ [ellipsoid,estr] = utmgeoid(utmstruct.zone);
+ %ellipsoid = almanac('earth','wgs84','meters');
+ %ellipsoid = almanac('earth','clarke66','meters');
+ %ellipsoid = [6.378206400000000e+06 0.082271854223002];
-r_dtor=1.74532925199d-2;
-i_ft=0;
-a_griddes=['C','D','E','F','G','H','J',...
-'K','L','M','N','P','Q','R','S','T','U',...
-'V','W','X'];
-r_a=6378206.4d0;
-r_e2=0.00676865799761d0;
-r_k0=0.9996d0; %scale at center
-r_lat0=0.0d0;
-r_fe=5e5;
-r_fn(1)=0;
-r_fn(2)=1e7; % N-S hemispheres??
+ ellipsoid = ellipsoid(1,:);
+ estr0 = strtrim(estr(1,:));
+ disp(sprintf(' ellipsoid %s : %.10e %.10e',estr0,ellipsoid));
+end
+
+utmstruct.geoid = ellipsoid; % assign ellipsoid
+utmstruct = defaultm(utmstruct);
+if i_type == 1
+ [y,x] = minvtran(utmstruct,xi,yi); % utm2ll
+else
+ [x,y] = mfwdtran(utmstruct,yi,xi); % ll2utm
+end
-% PROCESSING STEPS
+%------------------------------
-r_ep2 = r_e2/(1.d0 - r_e2);
-r_e4 = r_e2^2;
-r_e6 = r_e2^3;
-r_dtor = pi/180;
-
-if i_type == 0 % convert lat,lon to UTM
-
- xi=xi*r_dtor;
- yi=yi*r_dtor;
- if (i_zone==0)
- i_zone = fix(mod(xi+3.d0*pi,2.d0*pi)/(r_dtor*6.d0)) + 1;
- i_zone = max(min(i_zone,60),1);
- end
- r_lon0 = -pi + 6.d0*r_dtor*(i_zone-1) + 3.d0*r_dtor; % central meridian
-
- r_n = r_a/sqrt(1.d0 - r_e2*sin(yi)^2);
- r_t = tan(yi)^2;
- r_t2 = r_t^2;
- r_c = r_ep2*cos(yi)^2;
- r_ba = (xi - r_lon0)*cos(yi);
- r_a2 = r_ba^2;
- r_a3 = r_ba*r_a2;
- r_a4 = r_ba*r_a3;
- r_a5 = r_ba*r_a4;
- r_a6 = r_ba*r_a5;
- r_m = r_a*((1.0d0-r_e2/4 - 3.0d0*r_e4/64.0d0-5.d0*r_e6/256.d0)*yi - (3.d0*r_e2/8.d0 + 3.d0*r_e4/32.d0 +45.d0*r_e6/1024.d0)*sin(2.d0*yi) + (15.d0*r_e4/256.d0 +45.d0*r_e6/1024.d0)*sin(4.d0*yi) - (35.d0*r_e6/3072.d0)*sin(6.d0*yi));
- r_m0 = r_a*((1.d0-r_e2/4 - 3.d0*r_e4/64.d0 -5.d0*r_e6/256.d0)*r_lat0 - (3.d0*r_e2/8.d0 + 3.d0*r_e4/32.d0 +45.d0*r_e6/1024.d0)*sin(2.d0*r_lat0) + (15.d0*r_e4/256.d0 +45.d0*r_e6/1024.d0)*sin(4.d0*r_lat0) - (35.d0*r_e6/3072.d0)*sin(6.d0*r_lat0));
-
- r_v(1) = r_k0*r_n*(r_ba+(1.d0-r_t+r_c)*r_a3/6.d0 +(5.d0-18.d0*r_t+r_t2+72.d0*r_c-58.d0*r_ep2)*r_a5/120.d0);
- r_v(1) = r_v(1) + r_fe;
-
- r_v(2) = r_k0*(r_m - r_m0 + r_n*tan(yi)*( r_a2/2.d0 + (5.d0-r_t+9.d0*r_c+4.d0*r_c^2)*(r_a4/24.d0) + (61.d0-58.d0*r_t+r_t2+600.d0*r_c-330.d0*r_ep2)*(r_a6/720.d0) ));
- if yi >= 0
- r_v(2) = r_v(2) + r_fn(1);
- else
- r_v(2) = r_v(2) + r_fn(2);
- end
-
- r_k = r_k0*(1.d0+(1.d0+r_ep2*cos(yi)^2)*(r_v(1)-r_fe)^2/(2.d0*(r_k0^2)*r_n^2));
-
- i_gi = fix((yi/r_dtor+80.d0)/8.d0) + 1;
- i_gi = max(min(i_gi,20),1);
- a_grid = a_griddes(i_gi);
-
- xo=r_v(1);
- yo=r_v(2);
-
-elseif i_type == 1 % convert UTM to lat,lon
-
- r_vu(1) = xi - r_fe;
- r_vu(2) = yi;
- if r_vu(2) >= r_fn(2) % LOOK at this if in the S hemisphere!!
- r_vu(2) = yi - r_fn(2);
- end
- r_lon0 = -pi + 6.d0*r_dtor*(i_zone-1) + 3.d0*r_dtor;
-
- r_et = sqrt(1.d0-r_e2);
- r_e1 = (1.d0-r_et)/(1.d0+r_et);
- r_e12 = r_e1^2;
- r_e13 = r_e1*r_e12;
- r_e14 = r_e1*r_e13;
- r_m = r_vu(2)/r_k0;
- r_mu = r_m/(r_a*(1.d0-r_e2/4.d0-3.d0*r_e4/64.d0-5.d0*r_e6/256.d0));
- r_lat1 = r_mu + (3.d0*r_e1/2.d0-27.d0*r_e13/32.d0)*sin(2.d0*r_mu)+(21.d0*r_e12/16.d0-55.d0*r_e14/32.d0)*sin(4.d0*r_mu) +(51.d0*r_e13/96.d0)*sin(6.d0*r_mu) +(1097.d0*r_e14/512.d0)*sin(8.d0*r_mu);
-
- r_n = r_a/sqrt(1.d0 - r_e2*sin(r_lat1)^2);
- r_r = (r_a*(1.d0-r_e2))/sqrt(1.d0 - r_e2*sin(r_lat1)^2)^3;
- r_t = tan(r_lat1)^2;
- r_t2 = r_t^2;
- r_c = r_ep2*cos(r_lat1)^2;
- r_c2 = r_c^2;
- r_d = r_vu(1)/(r_n*r_k0);
- r_d2 = r_d^2;
- r_d3 = r_d2*r_d;
- r_d4 = r_d3*r_d;
- r_d5 = r_d4*r_d;
- r_d6 = r_d5*r_d;
-
- yo = r_lat1 - (r_n*tan(r_lat1)/r_r)*(r_d2/2.d0+(5.d0+3.d0*r_t+10.d0*r_c-4.d0*r_c2-9.d0*r_ep2)*r_d4/24.d0 +(61.d0+90*r_t+298.d0*r_c+45.d0*r_t2-252.d0*r_ep2-3.d0*r_c2)*(r_d6/720.d0));
- xo = r_lon0 + (r_d - (1.d0+2.d0*r_t+r_c)*r_d3/6.d0 +(5.d0-2.d0*r_c+28.d0*r_t-3.d0*r_c2+8.d0*r_ep2+24.d0*r_t2)*(r_d5/120.d0))/cos(r_lat1);
- xo=xo/r_dtor;
- yo=yo/r_dtor;
-
+if 0==1
+ % to view UTM zones using Matlab GUI
+ figure;
+ axesm utm
+ axesmui
+ % THEN CLICK "ZONE"
+
+ % to determine which zone a point is in, and the lat-lon bounds of the zone
+ p1 = [32 -118];
+ zs = utmzone(p1);
+ [zonelats zonelons] = utmzone(zs);
+ disp(sprintf('UTM zone %s extends from %.1f to %.1f in longitude and %.1f to %.1f in latitude',...
+ zs,zonelons(1),zonelons(2),zonelats(1),zonelats(2)));
+
+ % EXAMPLE
+ format long, clc, clear, close all
+ x0 = -118; y0 = 32;
+ %x0 = 93.45; y0 = 7.65;
+ p1 = [y0 x0];
+ zs = utmzone(p1);
+ [x1,y1] = utm2ll(x0,y0,zs,0);
+ [x2,y2] = utm2ll(x1,y1,zs,1);
+ disp(sprintf('original lon-lat point: %f, %f',x0,y0));
+ disp(sprintf('recovered lon-lat point: %f, %f',x2,y2));
end
-
-%==========================================================================
\ No newline at end of file
+
+%==========================================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m 2011-06-13 21:59:19 UTC (rev 18609)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m 2011-06-13 22:18:30 UTC (rev 18610)
@@ -1,29 +0,0 @@
-%
-% function [x,y] = utm2llvec(xi,yi,i_zone,i_type)
-% Carl Tape, 10-Nov-2004
-% printed xxx
-%
-%
-% calls utm2ll.m
-% called by omsplotsC.m
-%
-
-function [x,y] = utm2llvec(xi,yi,i_zone,i_type)
-
-[a,b] = size(xi);
-num = a*b;
-xvec0 = reshape(xi,num,1);
-yvec0 = reshape(yi,num,1);
-xvec = zeros(num,1);
-yvec = zeros(num,1);
-
-for ii=1:length(xi)
- %if mod(length(xi),1000)==0
- % disp(sprintf('UTM to LATLON : %i out of %i',ii,length(xi)));
- %end
- [xvec(ii),yvec(ii)] = utm2ll(xvec0(ii),yvec0(ii),i_zone,i_type);
-end
-x = reshape(xvec,a,b);
-y = reshape(yvec,a,b);
-
-%==========================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m 2011-06-13 21:59:19 UTC (rev 18609)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m 2011-06-13 22:18:30 UTC (rev 18610)
@@ -32,14 +32,17 @@
iplot_vert = 0;
ihorz = 0;
-iplot_horz = 0;
+iplot_horz = 1;
iplot_horz_coverage = 0;
-iplot_horz_ker = 1;
+iplot_horz_ker = 0;
-iwrite = 0;
+iwrite = 1;
ifigure = 0;
imask = 1; % write out masked versions with NaN (=1) or not (=0)
+% UTM zone
+szone = '11S';
+
%-----------------------------------
% get a set of plotting gridpoints to extract a function value at a
@@ -47,14 +50,14 @@
% bounds for SPECFEM3D simulations
ax0 = [-121.6 -114.7 32.2 36.8];
-ax0_utm = axes_utm2ll(ax0,11,0);
+ax0_utm = axes_utm2ll(ax0,szone,0);
%ax0_utm = [ 0.3271 0.4698 3.6792 3.7914] * 1e6; % zoom
xmin = ax0_utm(1); xmax = ax0_utm(2);
ymin = ax0_utm(3); ymax = ax0_utm(4);
% outline of SPECFEM region
[xbox_utm,ybox_utm] = boxpts(ax0_utm,500);
-[xbox_ll,ybox_ll] = utm2llvec(xbox_utm,ybox_utm,11,1);
+[xbox_ll,ybox_ll] = utm2ll(xbox_utm,ybox_utm,szone,1);
figure; plot(xbox_utm, ybox_utm, '.');
@@ -63,8 +66,8 @@
if ihorz == 1
- irun = 2; % 1 for 0 km, 1 km, 2 km, etc
- % 2 for 0.5 km, 1.5 km, etc
+ irun = 2; % 1 for 1 km depth increment
+ % 2 for 0.5 km depth increment
dir1 = [pdir 'horz_' sprintf('%2.2i',irun) '/'];
switch irun
@@ -90,8 +93,8 @@
N = length(Xvec);
% convert back to lat-lon
- disp('Entering utm2llvec.m...');
- [Xvec_ll,Yvec_ll] = utm2llvec(Xvec,Yvec,11,1);
+ disp('Entering utm2ll.m...');
+ [Xvec_ll,Yvec_ll] = utm2ll(Xvec,Yvec,szone,1);
figure; plot(Xvec_ll, Yvec_ll, '.','markersize',2);
disp(sprintf('horizontal mesh is %i (x) by %i (y) for %i points total',numx,numy,N));
@@ -118,6 +121,7 @@
for ii=1:ndep
zdep = depvec(ii)
filename = sprintf('horz_xc_%3.3i',ii);
+ disp(['writing file: ' dir1 filename]);
fid = fopen([dir1 filename],'w');
for kk=1:N
fprintf(fid,'%16.6e%16.6e%16.6e\n',Xvec(kk),Yvec(kk),zdep);
@@ -137,7 +141,7 @@
[rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(recfile);
stanet = merge_arrays(stnm,netwk,'.');
srcfile = '/home/carltape/results/SOURCES/socal_16/EIDs_lonlat_otime';
- [junk1,eid,junk2,junk3,sdep,junk5,junk6,slon,slat] = textread(srcfile,'%f%f%s%s%f%f%f%f%f');
+ [~,eid,~,~,sdep,~,~,slon,slat] = textread(srcfile,'%f%f%s%s%f%f%f%f%f');
% irun
% 1 src-rec
@@ -404,7 +408,7 @@
preseg = preseg_vec(ii);
postseg = postseg_vec(ii);
[lon_gc1,lat_gc1,rdist,azis] = raypath(lon1,lat1,lon2,lat2,preseg,postseg,dx); % KEY COMMAND
- [lon_gc1_utm,lat_gc1_utm] = utm2llvec(lon_gc1,lat_gc1,11,0);
+ [lon_gc1_utm,lat_gc1_utm] = utm2ll(lon_gc1,lat_gc1,szone,0);
% exclude the points on the ray that are outside the mesh,
% then convert back to lat-lon
@@ -413,7 +417,7 @@
azis = azis(inds);
lon_gc1_utm = lon_gc1_utm(inds);
lat_gc1_utm = lat_gc1_utm(inds);
- [lon_gc1,lat_gc1] = utm2llvec(lon_gc1_utm,lat_gc1_utm,11,1);
+ [lon_gc1,lat_gc1] = utm2ll(lon_gc1_utm,lat_gc1_utm,szone,1);
hnum = length(lon_gc1);
% compute distance to each of the six faults
@@ -454,7 +458,7 @@
% % get points on the ray path
% %[lat_gc1,lon_gc1] = track1(lat1,lon1,az1,km2deg(rng1),[],'degrees',hnum);
% [lat_gc1,lon_gc1] = gcwaypts(lat1,lon1,lat2,lon2,hnum-1);
-% [lon_gc1_utm,lat_gc1_utm] = utm2llvec(lon_gc1,lat_gc1,11,0);
+% [lon_gc1_utm,lat_gc1_utm] = utm2ll(lon_gc1,lat_gc1,szone,0);
%
% % compute the distance in km from one point to the end --
% % this is used in the 2D plots
@@ -484,6 +488,7 @@
if iwrite == 1
% xyz files
filename = sprintf('vert_%2.2i_xc_%3.3i_input',irun,ii);
+ disp(['writing file: ' dir1 filename]);
fid = fopen([dir1 filename],'w');
for kk=1:hnum % loop over horizontal points in the ray path
for jj=1:znum
@@ -566,20 +571,20 @@
[lon1,lat1,d1,z1,lon2,lat2,d2,z2,az,iflip] = textread(file0,'%f%f%f%f%f%f%f%f%f%f');
[lon_gc1,lat_gc1,lon_gc1_utm,lat_gc1_utm,rdist] = textread(file1,'%f%f%f%f%f');
[utm_x,utm_y,z,d] = textread(file2,'%f%f%f%f');
- [junk1,junk2,junk3,vs1,cdist] = textread(file3,'%f%f%f%f%f');
- [junk1,junk2,junk3,vs2,cdist] = textread(file4,'%f%f%f%f%f');
- [junk1,junk2,junk3,vb1,cdist] = textread(file5,'%f%f%f%f%f');
- [junk1,junk2,junk3,vb2,cdist] = textread(file6,'%f%f%f%f%f');
+ [~,~,~,vs1,cdist] = textread(file3,'%f%f%f%f%f');
+ [~,~,~,vs2,cdist] = textread(file4,'%f%f%f%f%f');
+ [~,~,~,vb1,cdist] = textread(file5,'%f%f%f%f%f');
+ [~,~,~,vb2,cdist] = textread(file6,'%f%f%f%f%f');
n = length(vs1);
%lon1 = lon_gc1(1); lat1 = lat_gc1(1);
%lon2 = lon_gc1(end); lat2 = lat_gc1(end);
% turn negative points into NaN values
- ineg = find(vs1 < 0); vs1(ineg) = NaN;
- ineg = find(vs2 < 0); vs2(ineg) = NaN;
- ineg = find(vb1 < 0); vb1(ineg) = NaN;
- ineg = find(vb2 < 0); vb2(ineg) = NaN;
+ vs1(vs1 < 0) = NaN;
+ vs2(vs2 < 0) = NaN;
+ vb1(vb1 < 0) = NaN;
+ vb2(vb2 < 0) = NaN;
% compute Vp
vp1 = zeros(n,1); vp2 = zeros(n,1);
@@ -670,11 +675,11 @@
if iplot_horz==1
colors;
- irun = 1;
+ irun = 2; % index into set of horizontal cross sections
stirun = sprintf('%2.2i',irun);
dir1 = [pdir 'horz_' stirun '/'];
smodel1 = 'm00';
- smodel2 = 'm01'; % m16, m01, etc
+ smodel2 = 'm16'; % m16, m01, etc
modlab = 'vs';
smodeltag = [smodel2 '_' smodel1];
@@ -689,8 +694,10 @@
cnorm = 10^cfac;
% load files needed for horizontal cross-sections
- [inds,depvec] = textread([dir1 'horz_' stirun '_xc_cuts'],'%f%f');
- [xlon,ylat,xutm,yutm] = textread([dir1 'horz_' stirun '_xc_lonlat'],'%f%f%f%f');
+ fname = [dir1 'horz_xc_cuts']; disp(['file : ' fname]);
+ [inds,depvec] = textread(fname,'%f%f');
+ fname = [dir1 'horz_xc_lonlat']; disp(['file : ' fname]);
+ [xlon,ylat,xutm,yutm] = textread(fname,'%f%f%f%f');
% for some reason matlab interpolation needs to work with smaller numbers,
% so we subtract a reference value
@@ -709,7 +716,7 @@
%plot(xlon(indmask),ylat(indmask),'ro');
%pmin = 1; pmax = length(depvec);
- pmin = 3; pmax = pmin;
+ pmin = 21; pmax = pmin;
%pmin = 1; pmax = 12;
for ip = pmin:pmax
@@ -720,9 +727,11 @@
tag1 = [modlab '_' smodel1];
tag2 = [modlab '_' smodel2];
file1 = [dir1 tag1 '/horz_xc_' tag1 '_' stip '.gmt'];
- [junk1,junk2,junk3,cval1,cdist] = textread(file1,'%f%f%f%f%f');
+ disp(['file :' file1]);
+ [~,~,~,cval1,cdist] = textread(file1,'%f%f%f%f%f');
file2 = [dir1 tag2 '/horz_xc_' tag2 '_' stip '.gmt'];
- [junk1,junk2,junk3,cval2,cdist] = textread(file2,'%f%f%f%f%f');
+ disp(['file :' file2]);
+ [~,~,~,cval2,cdist] = textread(file2,'%f%f%f%f%f');
if imask==1
% load coverage mask
@@ -730,13 +739,13 @@
%stip0=stip;
%if ip==3, stip0 = '002'; end
filem = [dir1 tagm '/horz_xc_' tagm '_' stip '.gmt'];
- [junk1,junk2,junk3,cvalm,cdist] = textread(filem,'%f%f%f%f%f');
+ [~,~,~,cvalm,cdist] = textread(filem,'%f%f%f%f%f');
end
% remove points above the topography -- these are set as negative
% numbers in the interpolation algorithm
- ineg = find(cval1 < 0); cval1(ineg) = NaN;
- ineg = find(cval2 < 0); cval2(ineg) = NaN;
+ cval1(cval1 < 0) = NaN;
+ cval2(cval2 < 0) = NaN;
% compute the difference -- GMT does not seem to like operating on
% NaN vectors, so we have to do it here and then write it out
@@ -792,7 +801,9 @@
if iwrite == 1
% write out a file with five columns: lon, lat, m00 (m/s), m16 (m/s), ln(m16/m00)
- fid = fopen([dir1 'horz_xc_' modlab '_' smodeltag '_' stip '_mask' num2str(imask) '.dat'],'w');
+ fname = [dir1 'horz_xc_' modlab '_' smodeltag '_' stip '_mask' num2str(imask) '.dat'];
+ disp(['writing file : ' fname]);
+ fid = fopen(fname,'w');
for kk=1:length(xlon)
fprintf(fid,'%16.6e%16.6e%16.6e%16.6e%16.6e\n',...
xlon(kk),ylat(kk),cval1(kk),cval2(kk),cval3(kk));
@@ -850,7 +861,7 @@
% load output from cluster code
tag1 = [modlab '_' smodeltag];
file1 = [dir1 tag1 '/horz_xc_' tag1 '_' stip '.gmt'];
- [junk1,junk2,junk3,cval1,cdist] = textread(file1,'%f%f%f%f%f');
+ [~,~,~,cval1,cdist] = textread(file1,'%f%f%f%f%f');
% for some reason matlab interpolation needs to work with smaller numbers,
% so we subtract a reference value
@@ -860,7 +871,7 @@
yutm = yutm - yutm0;
% remove points above the topography
- ineg = find(cval1 < 0); cval1(ineg) = NaN;
+ cval1(cval1 < 0) = NaN;
% normalize by number of windows or number of events
cval1 = cval1 / Nfac;
@@ -946,6 +957,8 @@
ftags = {'Kmu_rayleigh','Kmu_reflection_one','Kmu_reflection_two'};
% load files needed for horizontal cross-sections
+ fname = [dir1 'horz_' stirun '_xc_cuts'];
+ disp(['file : ' fname]);
[inds,depvec] = textread([dir1 'horz_' stirun '_xc_cuts'],'%f%f');
[xlon,ylat,xutm,yutm] = textread([dir1 'horz_' stirun '_xc_lonlat'],'%f%f%f%f');
@@ -962,11 +975,11 @@
% load output from cluster code
file1 = [dir1 'kernels/horz_xc_' ftags{1} '_' fsmooth '_' stip '.gmt'];
- [junk1,junk2,junk3,K1,junk4] = textread(file1,'%f%f%f%f%f');
+ [~,~,~,K1,~] = textread(file1,'%f%f%f%f%f');
file2 = [dir1 'kernels/horz_xc_' ftags{2} '_' fsmooth '_' stip '.gmt'];
- [junk1,junk2,junk3,K2,junk4] = textread(file2,'%f%f%f%f%f');
+ [~,~,~,K2,~] = textread(file2,'%f%f%f%f%f');
file3 = [dir1 'kernels/horz_xc_' ftags{3} '_' fsmooth '_' stip '.gmt'];
- [junk1,junk2,junk3,K3,junk4] = textread(file3,'%f%f%f%f%f');
+ [~,~,~,K3,~] = textread(file3,'%f%f%f%f%f');
% for some reason matlab interpolation needs to work with smaller numbers,
% so we subtract a reference value
More information about the CIG-COMMITS
mailing list