Search This Blog

Friday, July 6, 2012

Cruise specific sea surface temperature/salinity animations

I'm currently sitting offshore south of Iceland about to pass over the Mid-Atlantic Ridge.  We're shooting EM302 multibeam and, as always, I'm curious about what to expect from the ocean's temperature and salinity.

The RTOFS grids that I've been playing with for the past while have come in handy again.  I wrote a script to pull down the imagery that's generated on a daily basis on the RTOFS website and a bit of sorcery with ImageMagick's "convert" program massages it into a little animation for me to see the sea surface temperature and salinity every day.

Here's an example:


And here's the script:




#!/bin/bash


dir=`date '+%Y%m%d'`


rm -rf $dir
mkdir $dir
cd $dir


for field in temperature salinity; do


        echo "Doing field" $field


        for hour in 024 048; do
                echo "Doing nowcast" $hour
                wget http://polar.ncep.noaa.gov/global/nctest/images/large/rtofs_arctic_$field\_n$hour\_0.png
        done


        for hour in 024 048 072 096 120 144; do
                echo "Doing forecast" $hour
                wget http://polar.ncep.noaa.gov/global/nctest/images/large/rtofs_arctic_$field\_f$hour\_0.png
        done
done


for file in *.png; do


        prefix=`basename $file .png`
        convert -crop 1840x120+0+1131 $file scale.ppm


        # I have to convert to .ppm to toss out rotation info in png header such that 
        # the crop operation works?  Otherwise, the crop geometry is in the original unrotated
        # image geometry.
        convert -rotate 135 $file temp.ppm
        convert -crop 450x360+600+1000 temp.ppm temp2.ppm


        label=`echo $file | awk -F_ '{print '$dir',$3,$4}' `


        convert -draw "image Over 0,310 450,52 scale.ppm" -fill white -stroke black -pointsize 24 -draw "text 10,30 '$label'" temp2.ppm $file


done

for temperature_file in *temperature*png; do
        output_file=`echo $temperature_file | sed -e 's/temperature/TS/' `
        salinity_file=`echo $temperature_file | sed -e 's/temperature/salinity/' `

        convert -extent 900x360 -draw "image Over 450,0 450,360 $salinity_file" $temperature_file $output_file
done

convert -loop 0 -delay 50 rtofs_arctic_temperature_f*png temperature_animation.gif
convert -loop 0 -delay 50 rtofs_arctic_salinity_f*png salinity_animation.gif
convert -loop 0 -delay 50 rtofs_arctic_TS_f*png TS_animation.gif

rm temp.ppm temp2.ppm scale.ppm



Thursday, July 5, 2012

Splitting up Kongsberg Watercolumn files

Sometimes you just forget to save Kongsberg multibeam water column data into a separate file.  This just happened to me on a cruise and I found that the >2GB file sizes made my 32-bit software puke.  Luckily, the file sizes didn't exceed 4GB so I decided to write a file splitter in python the pulls apart the original .all file and outputs a new .all file purged of water column datagrams AND a separate .wcd file.  Here's my first cut at it, it's a script that you feed a list of filenames and it creates a "split" subdirectory for each file and writes a split .all/.wcd combination into the split subdirectory, for example, the file:

    20120702/0003_20120702_122222_FK_EM710.all

will split into:

    20120702/split/0003_20120702_122222_FK_EM710.all    
    20120702/split/0003_20120702_122222_FK_EM710.wcd

Give it a try, I hope it doesn't nuke your data.


#!/usr/bin/env python2.6


import os
import struct
import time
import sys 


file_count=0
debug=False


dir="split"


for filename in sys.argv:


    file_count += 1


    if (file_count == 1): 
        # I'm too lazy to parse command line args so just skipping over the 
        # script name (which is arg zero in the list)
        continue


    file = open(filename, 'rb')
    filesize = os.path.getsize(filename)


    # What is the path to the input file without the filename?
    filepath=os.path.dirname(filename)
    fileprefix=os.path.basename(filename)
    if debug or True:
        print "Doing file",filename
        print "Got file path",filepath
        print "Got file basename",fileprefix


    # Join the file's directory path with the usual output subdirectory name
    outdir=os.path.join(filepath,dir)


    if not os.path.exists(outdir):
        os.makedirs(outdir)


    split_allname = os.path.join(outdir, fileprefix)
    split_wcdname = split_allname.replace(".all",".wcd")


    if debug or True:
        print filename, "will split into", split_allname, split_wcdname



    if not os.path.exists(split_allname):
        split_allfile = open(split_allname,"wb")
        split_wcdfile = open(split_wcdname,"wb")
    else:
        print "Skipping", filename, "since it's already split!"
        file.close()
        continue


    last_percent = 0
    while True:


        # Make sure we don't try to read beyond the EOF
        if (file.tell() + 6 > filesize):
            break


        line = file.read(6)


        header = struct.unpack("


        rawlength=line[0:3]
        length = header[0]
        stx = header[1]
        id = header[2]


        if (stx != 2):
            if debug:
                print 'STX not found, trying next datagram at position',file.tell()-5
            file.seek(-5,1)
            continue


        if debug:
            print 'STX found, going to try for ETX now'


        # Make sure we don't try to read beyond the EOF
        if (file.tell() + (length-5) > filesize):
            file.seek(-5,1)
            continue


        file.seek(length-5,1)


        # Make sure we don't try to read beyond the EOF
        if (file.tell() + 3 > filesize):
            break


        line = file.read(3)
        footer = struct.unpack("
        etx = footer[0]
        checksum = footer[1]

        if (etx != 3):
            if debug:
                print 'ETX not found, trying next datagram at position',file.tell()-(length+3)
            file.seek(-(length+3),1)
            continue

        # Rewind to very beginning of the datagram, including the length field
        file.seek(-(length+4),1)
        data = file.read(length+4)

        if debug:
            print "Got id", id, "and length", length

        if (id == 0x49 or id == 0x69 or id == 0x52 or id == 0x55):
            # Stuff for both files
            split_allfile.write(data)
            split_wcdfile.write(data)
        elif (id == 0x6B):
            # Just for the watercolumn file
            split_wcdfile.write(data)
        else:
            # Everything else goes into the raw file
            split_allfile.write(data)

        percent=int(100.0 * file.tell()/filesize)

        if (percent%5 == 0 and percent != last_percent):
            print percent, "% done, ALL:",split_allfile.tell()," WCD:",split_wcdfile.tell()
            last_percent = percent

        if file.tell() >= filesize:
            break

    file.close()
    split_allfile.close()
    split_wcdfile.close()

print 'All done!'