IMMA files and diagnostics for Bear.

(Omitted 1884 because there is an unresolved doubling of the obs).

by_voyage/Bear/20CR.extract.ship.R

@@ -0,0 +1,60 @@
+# Get pressures from 20CR, for a selected ship.
+
+ship<-'Bear'
+
+library(GSDF.TWCR)
+library(parallel)
+
+# Get the observations for this ship
+o<-read.fwf(file='imma.annual/Bear.1911.imma', widths=c(4,2,2,4,5,6,36,5,4,5,12,5))
+o$Longitude<-o$V6/100
+o$Latitude<-o$V5/100
+o$Pressure<-as.numeric(o$V8)/10
+o$AT<-as.numeric(o$V10)/10
+o$SST<-as.numeric(o$V12)/10
+#w<-seq(1,1000) # Testing
+#o<-o[seq(1,3167),] # Eliminate the last few which would need 1947 data
+
+# Get mean and spread from 3.2.1 for each ob.
+
+get.comparisons<-function(i) {
+
+  year<-o$V1[i]
+  month<-o$V2[i]
+  day<-o$V3[i]
+  hour<-as.integer(o$V4[i]/100)
+  t2m<-TWCR.get.slice.at.hour('air.2m',year,month,day,hour,
+                              version='3.2.1',opendap=F)
+  tt<-GSDF.interpolate.ll(t2m,o$Latitude[i],o$Longitude[i])  
+  t2m<-TWCR.get.slice.at.hour('air.2m',year,month,day,hour,
+                              type='spread',
+                              version='3.2.1',opendap=F)
+  tt.spread<-GSDF.interpolate.ll(t2m,o$Latitude[i],o$Longitude[i])  
+  old<-TWCR.get.slice.at.hour('prmsl',year,month,day,hour,
+                              version='3.2.1',opendap=F)
+  mean<-GSDF.interpolate.ll(old,o$Latitude[i],o$Longitude[i])
+  old<-TWCR.get.slice.at.hour('prmsl',year,month,day,hour,
+                              type='spread',
+                              version='3.2.1',opendap=F)
+  spread<-GSDF.interpolate.ll(old,o$Latitude[i],o$Longitude[i])
+  #if(i==10) break
+  return(c(i,tt,tt.spread,mean,spread))
+}
+
+#lapply(seq_along(o$V1),get.comparisons)
+r<-mclapply(seq_along(o$V1),get.comparisons,mc.cores=6)
+r<-as.numeric(unlist(r))
+odr<-r[seq(1,length(r),5)]
+tt<-r[seq(2,length(r),5)]
+tt.spread<-r[seq(3,length(r),5)]
+mean<-r[seq(4,length(r),5)]
+spread<-r[seq(5,length(r),5)]
+
+# Output the result
+fileConn<-file(sprintf("obs.%s",ship))
+writeLines(sprintf("%d %d %d %d %f %f %f %f %f %f",
+                   o$V1[odr],o$V2[odr],o$V3[odr],o$V4[odr],
+                   o$Pressure[odr],mean/100,spread/100,
+                   o$AT[odr],tt-273.15,tt.spread),
+                   fileConn)
+close(fileConn)

by_voyage/Bear/O_R_plot.R

@@ -0,0 +1,90 @@
+# Plot a selected WD ship's pressure observations along with the
+#  reanalysis mean and spread along her route.
+
+ship<-'Bear'
+
+library(grid)
+library(chron)
+
+o<-read.table(sprintf("obs.%s",ship))
+o$Date<-chron(dates=sprintf("%04d/%02d/%02d",as.integer(o$V1),
+                                             as.integer(o$V2),
+                                             as.integer(o$V3)),
+             times=sprintf("%02d:30:30",as.integer(o$V4/100)),
+             format=c(dates = "y/m/d", times = "h:m:s"))
+tics<-pretty(o$Date)
+
+pdf(file=sprintf("%s_comparison.pdf",ship),
+    width=10,height=6,family='Helvetica',
+    paper='special',pointsize=12)
+
+pushViewport(viewport(width=0.98,height=0.5,x=0.01,y=0.01,
+                      just=c("left","bottom"),name="Page",clip='off'))
+   pushViewport(plotViewport(margins=c(4,4.5,1,0)))
+      pushViewport(dataViewport(o$Date,c(950,1050)))
+
+
+         grid.xaxis(at=as.numeric(tics),label=attr(tics,'label'),main=T)
+         grid.text('Date',y=unit(-3,'lines'))
+         grid.yaxis(main=T)
+         grid.text('Sea-level pressure (hPa)',x=unit(-4,'lines'),rot=90)
+
+
+         # Analysis spreads
+         gp=gpar(col=rgb(0.8,0.8,1,1),fill=rgb(0.8,0.8,1,1))
+         for(i in seq_along(o$V1)) {
+            x<-c(o$Date[i]-1/48,o$Date[i]+1/48,
+                 o$Date[i]+1/48,o$Date[i]-1/48)
+            y<-c(o$V6[i]-(o$V7[i])*2,
+                 o$V6[i]-(o$V7[i])*2,
+                 o$V6[i]+(o$V7[i])*2,
+                 o$V6[i]+(o$V7[i])*2)
+            grid.polygon(x=unit(x,'native'),
+                         y=unit(y,'native'),
+                      gp=gp)
+          }
+
+        # Observation
+         gp=gpar(col=rgb(0,0,0,1),fill=rgb(0,0,0,1))
+         grid.points(x=unit(o$Date,'native'),
+                     y=unit(o$V5,'native'),
+                     size=unit(0.005,'npc'),
+                     pch=20,
+                     gp=gp)
+      popViewport()
+   popViewport()
+upViewport()
+
+pushViewport(viewport(width=0.98,height=0.5,x=0.01,y=0.51,
+                      just=c("left","bottom"),name="Page",clip='off'))
+   pushViewport(plotViewport(margins=c(0,4.5,1,0)))
+      pushViewport(dataViewport(o$Date,c(-50,30)))
+
+
+         grid.yaxis(main=T)
+         grid.text('Air Temperature (C)',x=unit(-4,'lines'),rot=90)
+
+         # Analysis
+         gp=gpar(col=rgb(0.8,0.8,1,1),fill=rgb(0.8,0.8,1,1))
+         for(i in seq_along(o$V1)) {
+            x<-c(o$Date[i]-1/48,o$Date[i]+1/48,
+                 o$Date[i]+1/48,o$Date[i]-1/48)
+            y<-c(o$V9[i]-(o$V10[i])*2,
+                 o$V9[i]-(o$V10[i])*2,
+                 o$V9[i]+(o$V10[i])*2,
+                 o$V9[i]+(o$V10[i])*2)
+            grid.polygon(x=unit(x,'native'),
+                         y=unit(y,'native'),
+                      gp=gp)
+          }
+
+        # Observation
+         gp=gpar(col=rgb(0,0,0,1),fill=rgb(0,0,0,1))
+         grid.points(x=unit(o$Date,'native'),
+                     y=unit(o$V8,'native'),
+                     size=unit(0.005,'npc'),
+                     pch=20,
+                     gp=gp)
+      popViewport()
+   popViewport()
+popViewport()

by_voyage/Bear/clean_positions.perl

@@ -0,0 +1,54 @@
+#!/usr/bin/perl
+
+# Clean up the raw positions file for Bear.
+
+use strict;
+use warnings;
+
+#open(DIN,"<positions.out") or die;
+my @Lines = <STDIN>;
+#close(DIN);
+
+# Find all the dates for which there is a lat:lon position
+my %Ll;
+for(my $i=0;$i<scalar(@Lines);$i++) {
+    my $Line = $Lines[$i];
+    chomp($Line);
+    my @Fields = split /\t/,$Line;
+    if(!defined($Fields[0]) || $Fields[0] =~ /NA/) { next; }
+    if((defined($Fields[1]) && $Fields[1] =~ /(\d+)\D+(\d+)/) ||
+       (defined($Fields[2]) && $Fields[2] =~ /(\d+)\D+(\d+)/) ) {
+	$Ll{$Fields[0]}=1;
+    }
+}
+
+# Weed out duplicate dates for which there is a good Ll position, and
+#  infered port postitions which are not credible
+#open(DOUT,">positions.qc.out") or die;
+for(my $i=0;$i<scalar(@Lines);$i++) {
+    my $Line = $Lines[$i];
+    chomp($Line);
+    my @Fields = split /\t/,$Line;
+    if(defined($Ll{$Fields[0]})) {
+       if((defined($Fields[1]) && $Fields[1] =~ /(\d+)\D+(\d+)/) ||
+          (defined($Fields[2]) && $Fields[2] =~ /(\d+)\D+(\d+)/) ) {
+	   printf "%s\t%s\t%s\t      NA\t      NA\t      NA\n",
+               $Fields[0],$Fields[1],$Fields[2];
+       }
+       next;
+    }
+    else {
+     if($Fields[6] eq '                      tic' ||
+        $Fields[6] eq '                    horta' ||
+        $Fields[6] eq '                  behring' ||
+        $Fields[6] eq '                      goa' ||
+        $Fields[6] eq '                    horta' ||
+        $Fields[6] eq '                     nice' ||
+        $Fields[6] eq '                      gib' ||
+        $Fields[6] eq '                port said' ||
+        $Fields[6] eq '                   port t') { next; }
+     print "$Line\n";
+    }
+}  
+#close(DOUT);
+

by_voyage/Bear/makeall.sh

@@ -0,0 +1,18 @@
+# Make the IMMA file and diagnostic plots for Bear
+
+../../scripts/utilities/voyage_positions.perl --ship='Bear' > positions.out
+clean_positions.perl < positions.out > positions.weeded
+cd positions.qc.annual
+../split_by_year.perl < ../positions.weeded
+# Manually edited - remove dud positions from annual files
+cd ..
+../../scripts/utilities/voyage_t+p.perl --ship='Bear' > obs.out
+cd obs.qc.annual
+../split_by_year.perl < obs.out
+# Manually edited - remove dud positions from annual files
+# Manually edited - fix dud dates to obs.qc.out
+./to_imma.perl
+R --no-save < plot_voyage.R
+R --no-save < 20CR.extract.ship.R
+R --no-save < O_R_plot.R
+cp imma.annual/* ../../imma