timeserenu_nord.pl | searchcode

/geofunc/perlscripts/timeserenu_nord.pl

https://gitlab.com/gpskings/gpslibrary
Perl | 348 lines | 186 code | 55 blank | 107 comment | 18 complexity | 2acd8b6ec450273f55057972eac6004d MD5 | raw file

#!/usr/bin/perl -w
#
#
#
#
#-------------------------------------------------------
#
use strict;
use PDL;
use PDL::NiceSlice;
#use PDL::Core ':Internal';
use PDL::Slatec qw( matinv polyfit polfit );
use PDL::LinearAlgebra qw(mnorm);
use PDL::MatrixOps qw( inv );
use PDL::Fit::Polynomial;
use File::Basename qw(dirname basename fileparse);
use Date::Pcalc qw(Days_in_Year);
use Geo::Functions qw(rad_deg deg_rad);
use Geo::Constants qw(PI);

require "geo.pm";

#defining an initialising parameters
my $scale = '18/8';
#my $outdir = "kara/psfiles/itrf05/";
#my $outdir = "nord/psfiles/";
my $outdir = "/home/bgo/verkefni/gps_timar/disp/nord/psfiles/";

my ($paxis,$saxis, $area);
my ($year,$month,$day, $doy,$long,$lat,$height,$station,$N,$sigmamed,$ressq,$oversigma,$sc,$unctemp,$ind);
my ($filbase, $path,$suff,$psfileb,$psend,$psfile ); 
my ($start, $stop,$pdata,$coeff, $coeffuns, $label,$title,$nuvvel);
my (@start,@stop,@minmax,@minmaxt,@tmp,@int,@coeff,@coeffuns,@VDVA,@VDVB,@data,@xyz,@MEAN,@RMS);
my ($mean,$prms,$median,$min,$max,$adev,$rms,$days);

# scaling to mm or whatever
my $scaling = 1000;
my @comp = qw/East North Up/;
my @xlabel = ("","","Year");
my @pos = ("0.8 3.6", "5.6 0.6" );
my @sc = (5,4,3);


foreach  my $ARG (@ARGV) {
	 print $ARG,"\n";

	($filbase, $path,$suff ) = &fileparse($ARG,qw{enu.disp});
  $title = $filbase; #Station name
	#$psfileb = "test"; # psfile prefix 
	#$psend = "\.ps";
	#$psfile = ${psfileb}.${psend}; # psfile for gmt output
	$psfile = "${outdir}${filbase}enu\.ps"; # psfile for gmt output ####CHECK enu  and enu.disp make more general########3
	
	
	#changing variables through loop @tmp @minmax
 	
	#initiating plot
	system ("echo '0 0' | pstext  -JX${scale} -R-1/1.5/-1/3 -Y29.1 -P -K << EOF > ${psfile}\n0.3 -0.8 16 0 4 TC ${title}\nEOF\n");
	
	
	#reading in the data
	@data = rcols $ARG;
	
	#getting station information
	open FIL, $ARG;
	while (<FIL>) {
		 if (/^#{2}\s/) {
				chomp ((undef,undef,$xyz[0],$xyz[1],$xyz[2]) = split );
		 }
		 if (/^#{3}/) {
		 		chomp ((undef, $station, $long, $lat, $height ) = split);
				last;
		 }
		 
	}
	close FIL;

	
	my $interv = 0.01918;
	$start = &min($data[0]);
	$stop = &max($data[0]);
	$int[0] = which( $data[0] < 2007.0);
	$int[1] = setops( which( 2008.0 < $data[0]), 'AND', which( 2010.0 > $data[0]) ); 
	$int[2] = which( $data[0] > 2009.0);
	print $#int,"\n";
	#Fixing the Eurasian plate----------------------
	#Altamimi et al 2007 Absolute rotation pole fo Eurasia in ITRF2005 (56.330 deg,-95.979 deg,0.261 deg/my)
	#Altamimi et al 2007 Absolute rotation pole fo N. Amerika in ITRF2005 (-4.291 deg,-87.385 deg,0.192 deg/my)
	my @eurpole = map &rad_deg($_),(56.330,-95.979,0.261/1e6);
#	my @nAmerpole = map &rad_deg($_),(-4.291,-87.385,0.192/1e6);
	@eurpole = @{&ellecc([$eurpole[0],$eurpole[1]],[0,0,$eurpole[2]])};
#	@nAmerpole = @{&ellecc([$nAmerpole[0],$nAmerpole[1]],[0,0,$nAmerpole[2]])};
#	print join(" ",@eurpole),"\n";
#	print matinv(&xyz2enumatr(\@eurpole,1)) x ;
	$nuvvel = &eccell(&xyzell(&gdatum,\@xyz),&rotpole("",\@xyz,\@eurpole),1);
	print $nuvvel*1000;
#	my $namvel = &eccell(&xyzell(&gdatum,\@xyz),&rotpole("",\@xyz,\@nAmerpole),1);
#	print "-----------------\n";
	#$nuvvel = &eccell(&xyzell(&gdatum,\@xyz),&rotpole("EURA",\@xyz,),1);
	# my $vel=(pdl($namvel-$nuvvel)*$scaling)->squeeze ;
	#print $vel,"\n"; 
	#print mnorm($vel),"\n";
	#my $angle =  -(atan2($vel(0),$vel(1)))->sclr + &PI/2;
	#print  360+deg_rad($angle),"\n";
	#print "-----------------\n";
	$data[1] .= $data[1] - $nuvvel(0,1)*($data[0]-$start);
	$data[3] .= $data[3] - $nuvvel(0,0)*($data[0]-$start);
	#-------------------------------------------------------------------	

	for (my $i=0;$i<=2;$i++) {

		 # subtracting the average of first few days"
    If plate is passed
        Calculate the nnr-nuvela velocity of a given plate in geocentric coordinates.
    If plate is not passed
        Calculate the velocity in a given point in geocentric coordinates xyz with for a given rotation pole wxwywz
    """
    if plate != None:
        pass
    else:
        if xyz != None and wxwywz != None:
 
		 $days=0.002739726*7;
		 $ind = which($data[0] < $data[0]->at(0)+$days);
		 ($mean,$prms,$median,$min,$max,$adev,$rms) = 
		 stats($data[2*$i+1]->index($ind),1/($data[2*$i+2]->index($ind))**2);
		 $data[2*$i+1]=$data[2*$i+1]-$mean;

		 #the scale of axis
		 $paxis = "pa1.0f0.1:\"$xlabel[$i]\":/a20f5:\"$comp[$i]\"::.\"\":WSen"; #labels for axis
		 #		$saxis = "sa1of0.1"; #labels for axis

		 #data range
		 @minmaxt = ($data[0]->at(0)-0.2,$data[0]->at(-1)+0.2);
		 @minmax = (min($data[2*$i+1])*1000-5,max($data[2*$i+1])*1000+10);
		 $area = "$minmaxt[0]\/$minmaxt[1]\/$minmax[0]\/$minmax[1]"; #range
		 #----------------------------------------------------------------------------------------------	

		 #scaling of unsertainties----------------------------
		 #Fitting 
		 for (my $j =0 ; $j <= $#int ; $j++) { #fitting each year separetly
				#$unctemp = $data[2*$i+2];
				if ($j == 0) {
					 $data[2*$i+2] = $data[2*$i+2] * $sc[$i];
				}
				($coeff[$j], $coeffuns[$j]) = &smplsq($data[0]->index($int[$j]),
					 $data[2*$i+1]->index($int[$j]),
					 $data[2*$i+2]->index($int[$j]));
				($coeff[$j], $coeffuns[$j]) = ($coeff[$j]*${scaling},
					 $coeffuns[$j]*${scaling}); #changing to mm
				push @VDVB, ($coeff[$j],$coeffuns[$j]);

				#scaling of unsertainties----------------------------
				#if ($j == 0) {
				#($N) = dims($data[2*$i+1]->index(pdl[1..50]));
				#$sigmamed = sumover($data[2*$i+2]->index(pdl[1..50])*${scaling})/$N;
				#$ressq = ($data[2*$i+1]->index(pdl[1..50])*${scaling}-($data[0]->index(pdl[1..50])*($coeff[$j]->slice(1))+$coeff[$j]->slice(0)) )**2;
				#$oversigma = 1/(${scaling}*$data[2*$i+2]->index(pdl[1..50]));
#					print $sigmamed, " ", $ressq," ",$oversigma,"\n";

				#$sc = sqrt($N*sumover($ressq*$oversigma)/(($N-1)*sumover($oversigma)))/sqrt($sigmamed);
				#$unctemp = $data[2*$i+2]; 
				#$data[2*$i+2] = $unctemp * $sc;
				#print $sc,"\n";
#					print sqrt($ressq),"\n";
#					print $data[2*$i+2]*${scaling},"\n";
#					}

				#fiting again-----------------------------------------------------------------
#					($coeff[$j], $coeffuns[$j]) = &smplsq($data[0]->index($int[$j]),
#																					$data[2*$i+1]->index($int[$j]),
#																					$data[2*$i+2]->index($int[$j]));
#				  ($coeff[$j], $coeffuns[$j]) = ($coeff[$j]*${scaling},
#																		 $coeffuns[$j]*${scaling});
		 }	

		 #fitting everything

     my $cut = intersect(which($data[0] > 2007.1), which( $data[0] < 2007.7));
		 my $use = setops( which( 2010.3 > $data[0]), 'XOR' , $cut);
		 #my $use = which( 2010.3 > $data[0]);
		 ($coeff, $coeffuns) = &smplsq($data[0]->index($use),
				$data[2*$i+1]->index($use),
				$data[2*$i+2]->index($use));
		 ($coeff, $coeffuns) = ($coeff*${scaling},
				$coeffuns*${scaling}); #changing to mm
		 push @VDVA, ($coeff,$coeffuns);

		 #saving the average values of the last campaign
		 $ind = which($data[0] > 2008.0);
		 ($mean,$prms,$median,$min,$max,$adev,$rms) = 
		 stats($data[2*$i+1]->index($ind),1/($data[2*$i+2]->index($ind))**2);
		 push @MEAN, ($mean*${scaling},$rms*${scaling});

		 # preparing for ploting
		 @tmp = ([list $data[0]], [list $data[2*$i+1]], [list $data[2*$i+2]]);
		 $pdata="";
		 for (my $i=0;$i<=$#{$tmp[0]};$i++) { #input for psxy
				$pdata  .= sprintf "%s %s %s\n",  ${$tmp[0]}[$i], 
				${$tmp[1]}[$i]*${scaling}, 
				${$tmp[2]}[$i]*${scaling};
		 }

		 #ploting time series for each component
		 system ("psxy -JX$scale -R$area -B$paxis -Sc0.05 -W0.1/0 -G60 -Ey0.25c/2 -Y-9.01 -P -O -K << END >> ${psfile}\n${pdata}\nEND\n") == 0 or die "system failed: $?";

		 #plot the zero line
		 system("psxy -JX${scale} -R${area} -W2/0 -P -O -K << END >> ${psfile}\n$minmaxt[0] 0\n$minmaxt[1] 0\nEND\n");

		 #For plotting the fit	
		 for (my $j=0;$j<= $#int ;$j++) {
				@start = ($start-0.1, $coeff[$j]->at(0,0)+($start-0.1)*$coeff[$j]->at(1,0) );
				@stop = ($stop+0.1, $coeff[$j]->at(0,0)+($stop+0.1)*$coeff[$j]->at(1,0) );
				$label = sprintf( "%2.1f\261%2.1f mm/yr",$coeff[$j]->at(1,0),$coeffuns[$j]->at(1,0));

				#plot the fitt lines
				#if ($j == 0 ) {
				#					system("psxy -JX${scale} -R${area} -W2/0 -P -O -K << END >> ${psfile}\n$start[0] $start[1]\n$stop[0] $stop[1]\nEND\n");
				#put the fitt parameters on labels
				#			 system("pstext -JX${scale} -R0/8/0/4 -P -O -K << END >> ${psfile}\n$pos[$j] 13 0 0 5  $label\nEND\n") == 0 or die "system failed: $?";
				#					}
		 }
	}
	#finish the plot
	system "echo '0 0' | psxy -JX -R -Ss0.1 -G255 -P -O >> ${psfile}";

	my $file1 = "NORD0509enu.vel";
	#open for appending
	open FIL1, '>>' , $file1;
	# the fit for the hole period
	printf FIL1 "%2.6f %2.6f\t%2.4f %2.4f\t%2.4f %2.4f\t%2.4f %2.4f\t%s\n",
	$long,$lat,
	$VDVA[0]->at(1,0),$VDVA[1]->at(1,0),
	$VDVA[2]->at(1,0),$VDVA[3]->at(1,0),
	$VDVA[4]->at(1,0),$VDVA[5]->at(1,0),
	$station;
	#open for appending
	close FIL1;


	my $file2 = "NORD0506enu.vel";
	#open for appending
	open FIL2, '>>' , $file2;
	#printing the first part of the period
	 printf FIL2 "%2.6f %2.6f\t%2.4f %2.4f\t%2.4f %2.4f\t%2.4f %2.4f\t%s\n",
																						 							$long,$lat,
																 $VDVB[0]->at(1,0),$VDVB[1]->at(1,0),
																 $VDVB[6]->at(1,0),$VDVB[7]->at(1,0),
																 $VDVB[12]->at(1,0),$VDVB[13]->at(1,0),
																														$station;
	close FIL2;

	my $file3 = "NORD0809enu.vel";
	#open for appending
	open FIL3, '>>' , $file3;
	#printing the last part of the period
	printf FIL3 "%2.6f %2.6f\t%2.4f %2.4f\t%2.4f %2.4f\t%2.4f %2.4f\t%s\n",
																													$long,$lat,
																 $VDVB[2]->at(1,0),$VDVB[3]->at(1,0),
																 $VDVB[8]->at(1,0),$VDVB[9]->at(1,0),
															 $VDVB[14]->at(1,0),$VDVB[15]->at(1,0),
																														$station;

	my $file4 = "NORD0910enu.vel";
	#open for appending
	open FIL4, '>>' , $file4;
	#printing the last part of the period
	printf FIL4 "%2.6f %2.6f\t%2.4f %2.4f\t%2.4f %2.4f\t%2.4f %2.4f\t%s\n",
																													$long,$lat,
																 $VDVB[4]->at(1,0),$VDVB[5]->at(1,0),
																 $VDVB[10]->at(1,0),$VDVB[11]->at(1,0),
															 $VDVB[16]->at(1,0),$VDVB[17]->at(1,0),
																														$station;


	$year = 2007;
	$doy = 243;
	#	$day = $year+$doy/Days_in_Year($year,12); #calculating the difference at day $day
	$day = $stop; #calculating the difference at day $day
	my $Dday = ($day - $start)/2; 
	
	
	print join(" ",@MEAN),"\n";
	print $MEAN[4]-($VDVB[8]->(0)+$VDVB[8]->(1)*$day),"\n";

	$file4 = "NORDenu.dsp";
	#open for appending
	open FIL4, '>>' , $file4;
	#printing the last part of the period
						  printf FIL4 "%2.6f %2.6f\t%2.4f %2.4f\t%2.4f %2.4f\t%2.4f %2.4f\t%s\n",
																																		 $long,$lat,
																															$MEAN[0],$MEAN[1],
																															$MEAN[2],$MEAN[3],
	 														list($MEAN[4]-($VDVB[8]->(0)+$VDVB[8]->(1)*$day)),
										 list( sqrt(($MEAN[5] )**2 +  ( $VDVB[9]->(1)*$Dday )**2 )),
				 																															 $station;

#						  printf FIL4 "%2.6f %2.6f\t%2.4f %2.4f\t%2.4f %2.4f\t%2.4f %2.4f\t%s\n",
#																																		 $long,$lat,
#		list(($VDVB[2]->(0)+$VDVB[2]->(1)*$day)-($VDVB[0]->(0)+$VDVB[0]->(1)*$day)),
#					  list( sqrt( ($VDVB[3]->(1)*$Dday)**2 + ( $VDVB[1]->(1)*$Dday)**2 )),
#	 list( ($VDVB[6]->(0)+$VDVB[6]->(1)*$day)-($VDVB[4]->(0)+$VDVB[4]->(1)*$day)),
#						list(sqrt( ( $VDVB[7]->(1)*$Dday)**2 + ( $VDVB[5]->(1)*$Dday)**2 )),
#	list(($VDVB[10]->(0)+$VDVB[10]->(1)*$day)-($VDVB[8]->(0)+$VDVB[8]->(1)*$day)),
#				 list( sqrt(($VDVB[11]->(1)*$Dday )**2 +  ( $VDVB[9]->(1)*$Dday )**2 )),
#				 																															 $station;
	close FIL4;
	@MEAN = ();
  @VDVB = ();
	@VDVA = ();
}







# subroutines ------------------------------------------------------------


#----------------------------------------------------------------
#
# Least squere fit
#
#
#
#------------------------------------------------------------------

sub smplsq { 
	 my $M = zeroes(2,2);
	 my $V = zeroes(1,2);
	 my $X = shift();
	 my $Y = shift();
	 my $chi = shift();

	 $M(0,0) .= sumover($chi**-2); 
	 $M(0,1) .= sumover($X*$chi**-2); 
	 $M(1,0) .= $M(0,1);
	 $M(1,1) .= sumover($X**2*$chi**-2);
	 $V(0,0) .= sumover($Y*$chi**-2);
	 $V(0,1) .= sumover($Y*$X*$chi**-2);
	 if ($M(0,0)*$M(1,1)-$M(0,1)*$M(1,0)) {
			return ( transpose(matinv($M) x $V), sqrt(matinv($M)->diagonal(0,1)) );
	 } else {
			return pdl(0 , 0), pdl(0,0);
	 }
}