root/codar2shape/trunk/codar2shape/codar2shape_nc_radial.py

#
# PROGRAM:      codar2shape_nc_radial.py
#
# USAGE:        python codar2shape_nc_radial.py site_location
#              (site_locations are duck or haty)
#               Note:  User should edit variables under the 'Default / user 
#                      defined variables' section below to set desired processing
#                      parameters.  Variables include:
#                        url - OPeNDAP server (dataset) URL
#                        workspace - filesystem directory location for processing
#                        ShapeFile - name of output shapefile
#                        subSample - data subsampling (skip) factor (0 = full dataset)
#                        xUnit - X dimension name (ie. 'lon' or longitude)
#                        yUnit - Y dimension name (ie. 'lat' or latitude)
#                        tUnit - Time dimension name (ie. 'time_series'), if present.
#                        uUnit - U dimension name (ie. 'sstb' or Sea Surface Temp.)
#                        vUnit - V dimension name
#                        tValue - Specify which timestep in the time array to download.
#                                 If time dimension is not present in the dataset, the
#                                 program will ignore this parameter.
#                                 (default = -9999; program will automatically calculate
#                                 and use the last timestep in the time array)
#                        
    
#
# Default / user defined variables
#

workspace = "/home/jcleary/data/codar/nccoos/"
subSample = 0
xUnit = "longitude"
yUnit = "latitude"
tUnit = "time"
tValue = -9999 # -9999 specifies to use last timestep in the time array (default action)
uUnit = "radial_water_u"
vUnit = "radial_water_v"
lUnit = "location"
timeDimPresent = 0
noData = '-9999.0'

#
# Import required modules
#

import sys, time, os, shutil, string

from math import atan2
from math import pi

from dap.client import *
from dap.proxy import *
from types import *
import dap
from dap.parsers import dds, das

#
# arguments for individual CODAR sites
#
try:
    arg = sys.argv[1]
    if arg == "duck":
        url = "http://nemo.isis.unc.edu/cgi-bin/nph-dods/data/nc-coos/latest_v2.0/nccoos-duck-codar-latest.nc"
        outPrefix = "duck"
    elif arg == "haty":
        url = "http://nemo.isis.unc.edu/cgi-bin/nph-dods/data/nc-coos/latest_v2.0/nccoos-haty-codar-latest.nc"
        outPrefix = "haty"
except:
    print "No site argument supplied (options are duck or haty)."
    sys.exit()

#
# set pathnames based on site name
asciiOut = workspace + "ascii_" + outPrefix + ".csv"
ShapeOut = workspace + outPrefix + ".shp"



#-----------------------------------------------------------------------------
#
# Functions
#

def cleanWorkspace():                 # Sync and close (and delete) data files
    print
    try:
      asciiExport.close()
    except:
      print "ASCII data file not closed"
    try:
      os.remove(asciiOut)
    except:
      print "Unable to remove ASCII out file"
    print "Workspace cleaned, all temporary files removed."

#-----------------------------------------------------------------------------
#
# Main Function
#

print
print "***"
print "*** Importing OPeNDAP Codar radial vectors (point locations) to a Shapefile"
print "***"
print
print "Importing OPeNDAP data from URL: " + url

#
# Assign/validate X, Y, T, U & V vector variables, set X, Y, & T array sizes

data = open(url)  #JC

try:
  xVar = data[xUnit]
  xShape = xVar.shape[0]
except:
  print "Error,",xUnit,"dimension does not exist in",url
  sys.exit()

try:
  yVar = data[yUnit]
  yShape = yVar.shape[0]
except:
  print "Error,",yUnit,"dimension does not exist in",url
  sys.exit()

try:
  tVar = data[tUnit] 
  tShape = tVar.shape[0]
  timeDimPresent = 1
except:
  print "Note,",tUnit,"dimension does not exist in",url

try:
  uVar = data[uUnit] 
except:
  print "Error,",uUnit,"dimension does not exist in",url
  sys.exit()

try:
  vVar = data[vUnit]
except:
  print "Error,",vUnit,"dimension does not exist in",url
  sys.exit()


print
print "Subsample factor:",subSample
print "Time dim. present:",timeDimPresent

#
# Check the validity of the L & T dimensions in the U attrib. vector
print "U attibute type:",uUnit
a = 0
while 1:
  try:
    uVar.dimensions[a]
  except:
    break
  print uUnit,"dimension",a,":",uVar.dimensions[a]
  if uVar.dimensions[a] not in [lUnit, tUnit]:
    print "Error - " + uVar.dimensions[a] + " is an invalid dimension."
    sys.exit()
  a += 1

# Check the validity of the L & T dimensions in the V attrib. vector
print "V attibute type:",vUnit
b = 0
while 1:
  try:
    vVar.dimensions[b]
  except:
    break
  print vUnit,"dimension",b,":",vVar.dimensions[b]
  if vVar.dimensions[b] not in [lUnit, tUnit]:
    print "Error - " + vVar.dimensions[b] + " is an invalid dimension."
    sys.exit()
  b += 1  
  

#
# Establish time step value and time array index location
if timeDimPresent:                  # Deal w/ time dimension if present
  tArray = tVar[0:tShape]
  if (tValue <= noData):         # No time step specified by user
    tStep = tShape - 1              # Select single element (last) in time step array
    tValue = tArray[tShape - 1]
  else:                             # Validate user timestep value & set array location
    tStep = 0
    for tElement in tArray:
      if tElement == tValue:
        break
      tStep += 1
    else:
      print "Error: ", tValue, " is an invalid timestep value."
      print "Valid timestep values: ", tArray
      sys.exit()
  print "X array shape = ", xShape
  print "Y array shape = ", yShape
  print "time array step = ", tStep
  print "time step value = ", tValue
  timeHandle = time.strftime('%m-%d-%Y %H:00 UTC',(time.gmtime(tValue))) 
  print "Readable time stamp: ", timeHandle
  

#
# Download X, Y, U, & V data arrays from OPeNDAP server
print
print "Downloading " + xUnit + ", " + yUnit + ", " + lUnit + ', ' + uUnit + ", " + vUnit + " data arrays..."
xArray = xVar[0:xShape-1:subSample]
yArray = yVar[0:yShape-1:subSample]

if timeDimPresent:                  # 3 dimensional vector (T,L) for U and V arrays
  uArray = uVar[tStep, 0:xShape-1:subSample]
  vArray = vVar[tStep, 0:xShape-1:subSample]
  
else:                               # 2 dimensional vector (X,Y) for U and V arrays
  uArray = uVar[0:xShape-1:subSample, 0:yShape-1:subSample]
  vArray = vVar[0:xShape-1:subSample, 0:yShape-1:subSample]

print
print "Number of Long. values:",xShape
print "Number of Lat. values:",yShape
print "ASCII temp file:",asciiOut
print "Total number of dimensions in",uUnit,"=",a
print "Total number of dimensions in",vUnit,"=",b

#
# Setup temp. process file(s)
try:
  asciiExport = file(asciiOut,'w')
except IOError:
  print "Unable to create temp ASCII data file.  Check dir. rights..."
  sys.exit()


#
# Load data into ASCII data file (comma delimited -> ESRI shapefile - points)  
l = 0
t = 0
                                    # Write column header
asciiExport.write("time," + xUnit + "," + yUnit + "," + lUnit + ",U,V,speed_ms,speed_mph,speed_knots,dir,dir_ms\n")
                                    # Write X, Y, Z data to temp. ASCII file
print
print "Writing Timestamp, " + xUnit + ", " + yUnit + ", " + lUnit + ', ' + uUnit + ", " + vUnit + ", speed(3), and direction(2) data to " + asciiOut + "..."
while l < len(xArray):
    if timeDimPresent:              # 3 dimensional vector (T,X,Y) for U and V array
      asciiExport.write(timeHandle + ",")
      if float(xArray[l]) > 180.:   # shift X values from 0 - 360 to -180 - 180
        asciiExport.write(str(xArray[l] - 360.) + ",")
      else:
        asciiExport.write(str(xArray[l]) + ",")
      asciiExport.write(str(yArray[l]) + ",")
      asciiExport.write(str(l) + ",")
      asciiExport.write(str(uArray[t,l]) + ",")
      asciiExport.write(str(vArray[t,l]) + ",")

      # code to derive Speed from U and V components
      # (U^2 + V^2)^0.5
      speed = (uArray[t,l]**2 + vArray[t,l]**2)**.5  #in CM/S
      speed_ms = speed/100
      speed_mph = (speed/100) * 2.237
      speed_knots = (speed/100) * 1.944 
      asciiExport.write(str(speed_ms) + "," + str(speed_mph) + "," + str(speed_knots) + ",")

      # code to derive Direction from U and V components
      # atan2(-U,-V) * (180/pi) - 180
      direction = atan2(-uArray[t,l],-vArray[t,l]) * (180/pi) - 180
      if direction < 0:
        azimuth = direction + 360
      elif direction > 360:
        azimuth = direction - 360
      ms_azimuth = -azimuth 
      asciiExport.write(str(azimuth) + "," + str(ms_azimuth))
      asciiExport.write('\n')

    else:                           # 2 dimensional vector (X,Y) for U and V array
      if float(xArray[x]) > 180.:   # shift X values from 0 - 360 to -180 - 180
        asciiExport.write(str(xArray[l] - 360.) + ",")
      else:
        asciiExport.write(str(float(xArray[l])) + ",")
      asciiExport.write(str(yArray[l]) + ",")
      asciiExport.write(str(uArray[l]) + ",")
      asciiExport.write(str(vArray[l]))
      asciiExport.write('\n')
    l += 1

asciiExport.write(timeHandle + ",5,10,NaN,777,777,NaN,NaN,NaN,NaN,Nan")
asciiExport.close()

#
# OGR2OGR Process: ASCII CSV to Shapefile
print
print "Removing existing files (OGR cannot overwite shapefiles)"
print
os.system("rm -f " + workspace + outPrefix + "*")

print "Using OGR to create shapefile..."
cmd = '/usr/local/bin/ogr2ogr -f "ESRI Shapefile" -where \'U != "' + noData + '"\' ' + ShapeOut + ' /home/jcleary/data/source_files/' + outPrefix + '_pydap.vrt'
print cmd
os.system(cmd)

cleanWorkspace()

print
print "***"
print "*** Data import and conversion complete."
print "***"
print

sys.exit()