root/codar2shape/trunk/codar2shape/codar2shape_nc_totals.py

# --------------------------------------------------------------------------------
# PROGRAM:      codar2shape_nc_totals.py
#
# PURPOSE:      Python routine for direct import of OPeNDAP vector
#               datasets (mainly HF_RADAR) into CSV and then OGR2OGR into shapefile (vector) format.
#
# USAGE:        python codar2shape_nc_totals.py  
#               
#               Note:  User should edit variables under the 'Default / user 
#                      defined variables' section below to set desired processing
#                      parameters.  Variables include:
#                        workspace - filesystem directory location for processing
#                        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.
#                        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/seacoos/"
subSample = 0
xUnit = "lon"
yUnit = "lat"
tUnit = "time"
tValue = -9999 # -9999 specifies to use last timestep in the time array (default action)
timeDimPresent = 0 
arcNoData = -9999 
asciiOut = workspace + "ascii_totals.csv"
url = "http://nemo.isis.unc.edu/cgi-bin/nph-dods/data/nc-coos/latest_v2.0/nccoos-ouba-codar-latest.nc"
uUnit = "water_u"
vUnit = "water_v"
# add parsing of nodata value from DODS server to complete the -where clause
noData = "-9999.0"
outPrefix = "UNC"

#
# Import required modules
import sys, os, shutil, string, time
from time import time, mktime, gmtime, localtime, strftime
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

    
#------------------------------------------------------------------------------------------------------
#
# 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 data file"
    print "Workspace cleaned, all temporary files removed."

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

print
print "***"
print "*** Importing OPeNDAP HF Radar data streams to 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) 

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 X, Y, & 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 [xUnit, yUnit, tUnit]:
    print "Error - " + uVar.dimensions[a] + " is an invalid dimension."
    sys.exit()
  a += 1

# Check the validity of the X, Y, & 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 [xUnit, yUnit, 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 <= arcNoData):         # 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

#
# Use time(GMT) in shapefile file name
time_handle = strftime('%Y_%m_%d_%H_00_00',(gmtime(tValue))) 
formatTime = strftime('%m-%d-%Y %H:00 UTC',(gmtime(tValue)))
ShapeOut = workspace + "ouba.shp"

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

if timeDimPresent:                  # 3 dimensional vector (T,X,Y) for U and V arrays
  uArray = uVar[tStep, 0:xShape:subSample, 0:yShape:subSample]
  vArray = vVar[tStep, 0:xShape:subSample, 0:yShape:subSample]
  
else:                               # 2 dimensional vector (X,Y) for U and V arrays
  uArray = uVar[0:xShape:subSample, 0:yShape:subSample]
  vArray = vVar[0:xShape:subSample, 0:yShape: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 -> shapefile - points)  
x = 0
y = 0
t = 0
                                    # Write column header
asciiExport.write("inst,time," + xUnit + "," + yUnit + ",U,V,speed_ms,speed_mph,speed_knots," + "dir,dir_ms\n")
                                    # Write X, Y, Z data to temp. ASCII file
print
print "Writing institution, time, " + xUnit + ", " + yUnit + ", " + uUnit + ", " + vUnit + ", speed(3), and direction(2) data to " + asciiOut + "..."
while y < len(yArray):
  while x < len(xArray):
    if timeDimPresent:              # 3 dimensional vector (T,X,Y) for U and V array
      asciiExport.write(outPrefix + "," + time_handle + ",")
      if float(xArray[x]) > 180.:   # shift X values from 0 - 360 to -180 - 180
        asciiExport.write(str(xArray[x] - 360.) + ",")
      else:
        asciiExport.write(str(xArray[x]) + ",")
      asciiExport.write(str(yArray[y]) + ",")          
      asciiExport.write(str(uArray[t,x,y]) + ",")
      asciiExport.write(str(vArray[t,x,y]) + ",")
      #
      # code to derive Speed from U and V components
      # (U^2 + V^2)^0.5
      speed = (uArray[t,x,y]**2 + vArray[t,x,y]**2)**.5
      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,x,y],-vArray[t,x,y]) * (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
      asciiExport.write(outPrefix + "," + time_handle + ",")
      if float(xArray[x]) > 180.:   # shift X values from 0 - 360 to -180 - 180
        asciiExport.write(str(xArray[x] - 360.) + ",")
      else:
        asciiExport.write(str(float(xArray[x])) + ",")
      asciiExport.write(str(yArray[y]) + ",")
      asciiExport.write(str(uArray[x,y]) + ",")
      asciiExport.write(str(vArray[x,y]))
      asciiExport.write('\n')
    x += 1
  x = 0
  y += 1

asciiExport.write(outPrefix + "," + formatTime + ",5,10,777,777,NaN,NaN,NaN,NaN,NaN")  # bogus data record for TS labeling in MapServer
asciiExport.close()

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

#
# MV to NCCOOS folder
os.system("mv /home/jcleary/data/codar/seacoos/ouba.shp /home/jcleary/data/codar/nccoos/ouba.shp")   
os.system("mv /home/jcleary/data/codar/seacoos/ouba.dbf /home/jcleary/data/codar/nccoos/ouba.dbf")
os.system("mv /home/jcleary/data/codar/seacoos/ouba.shx /home/jcleary/data/codar/nccoos/ouba.shx")
os.system("mv /home/jcleary/data/codar/seacoos/ouba.prj /home/jcleary/data/codar/nccoos/ouba.prj")
print
print "Move UNC file to be used on NCCOOS.ORG"


cleanWorkspace()

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

sys.exit()