Index: codar2shape/trunk/codar2shape/codar2shape_nc_totals.py
===================================================================
---  (revision )
+++ codar2shape/trunk/codar2shape/codar2shape_nc_totals.py (revision 48)
@@ -1,0 +1,293 @@
+# --------------------------------------------------------------------------------
+# 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()