1 |
# -------------------------------------------------------------------------------- |
---|
2 |
# PROGRAM: codar2shape_nc_totals.py |
---|
3 |
# |
---|
4 |
# PURPOSE: Python routine for direct import of OPeNDAP vector |
---|
5 |
# datasets (mainly HF_RADAR) into CSV and then OGR2OGR into shapefile (vector) format. |
---|
6 |
# |
---|
7 |
# USAGE: python codar2shape_nc_totals.py |
---|
8 |
# |
---|
9 |
# Note: User should edit variables under the 'Default / user |
---|
10 |
# defined variables' section below to set desired processing |
---|
11 |
# parameters. Variables include: |
---|
12 |
# workspace - filesystem directory location for processing |
---|
13 |
# subSample - data subsampling (skip) factor (0 = full dataset) |
---|
14 |
# xUnit - X dimension name (ie. 'lon' or longitude) |
---|
15 |
# yUnit - Y dimension name (ie. 'lat' or latitude) |
---|
16 |
# tUnit - Time dimension name (ie. 'time_series'), if present. |
---|
17 |
# tValue - Specify which timestep in the time array to download. |
---|
18 |
# If time dimension is not present in the dataset, the |
---|
19 |
# program will ignore this parameter. |
---|
20 |
# (default = -9999; program will automatically calculate |
---|
21 |
# and use the last timestep in the time array) |
---|
22 |
|
---|
23 |
# |
---|
24 |
# Default / user defined variables |
---|
25 |
workspace = "/home/jcleary/data/codar/seacoos/" |
---|
26 |
subSample = 0 |
---|
27 |
xUnit = "lon" |
---|
28 |
yUnit = "lat" |
---|
29 |
tUnit = "time" |
---|
30 |
tValue = -9999 # -9999 specifies to use last timestep in the time array (default action) |
---|
31 |
timeDimPresent = 0 |
---|
32 |
arcNoData = -9999 |
---|
33 |
asciiOut = workspace + "ascii_totals.csv" |
---|
34 |
url = "http://nemo.isis.unc.edu/cgi-bin/nph-dods/data/nc-coos/latest_v2.0/nccoos-ouba-codar-latest.nc" |
---|
35 |
uUnit = "water_u" |
---|
36 |
vUnit = "water_v" |
---|
37 |
# add parsing of nodata value from DODS server to complete the -where clause |
---|
38 |
noData = "-9999.0" |
---|
39 |
outPrefix = "UNC" |
---|
40 |
|
---|
41 |
# |
---|
42 |
# Import required modules |
---|
43 |
import sys, os, shutil, string, time |
---|
44 |
from time import time, mktime, gmtime, localtime, strftime |
---|
45 |
from math import atan2 |
---|
46 |
from math import pi |
---|
47 |
from dap.client import * |
---|
48 |
from dap.proxy import * |
---|
49 |
from types import * |
---|
50 |
import dap |
---|
51 |
from dap.parsers import dds, das |
---|
52 |
|
---|
53 |
|
---|
54 |
#------------------------------------------------------------------------------------------------------ |
---|
55 |
# |
---|
56 |
# Functions |
---|
57 |
# |
---|
58 |
|
---|
59 |
def cleanWorkspace(): # Sync and close (and delete) data files |
---|
60 |
print |
---|
61 |
try: |
---|
62 |
asciiExport.close() |
---|
63 |
except: |
---|
64 |
print "ASCII data file not closed" |
---|
65 |
try: |
---|
66 |
os.remove(asciiOut) |
---|
67 |
except: |
---|
68 |
print "Unable to remove ASCII data file" |
---|
69 |
print "Workspace cleaned, all temporary files removed." |
---|
70 |
|
---|
71 |
#----------------------------------------------------------------------------- |
---|
72 |
# |
---|
73 |
# Main Function |
---|
74 |
# |
---|
75 |
|
---|
76 |
print |
---|
77 |
print "***" |
---|
78 |
print "*** Importing OPeNDAP HF Radar data streams to shapefile" |
---|
79 |
print "***" |
---|
80 |
print |
---|
81 |
print "Importing OPeNDAP data from URL: " + url |
---|
82 |
|
---|
83 |
# |
---|
84 |
# Assign/validate X, Y, T, U & V vector variables, set X, Y, & T array sizes |
---|
85 |
data = open(url) |
---|
86 |
|
---|
87 |
try: |
---|
88 |
xVar = data[xUnit] |
---|
89 |
xShape = xVar.shape[0] |
---|
90 |
except: |
---|
91 |
print "Error,",xUnit,"dimension does not exist in",url |
---|
92 |
sys.exit() |
---|
93 |
|
---|
94 |
try: |
---|
95 |
yVar = data[yUnit] |
---|
96 |
yShape = yVar.shape[0] |
---|
97 |
except: |
---|
98 |
print "Error,",yUnit,"dimension does not exist in",url |
---|
99 |
sys.exit() |
---|
100 |
|
---|
101 |
try: |
---|
102 |
tVar = data[tUnit] |
---|
103 |
tShape = tVar.shape[0] |
---|
104 |
timeDimPresent = 1 |
---|
105 |
except: |
---|
106 |
print "Note,",tUnit,"dimension does not exist in",url |
---|
107 |
|
---|
108 |
try: |
---|
109 |
uVar = data[uUnit] |
---|
110 |
except: |
---|
111 |
print "Error,",uUnit,"dimension does not exist in",url |
---|
112 |
sys.exit() |
---|
113 |
|
---|
114 |
try: |
---|
115 |
vVar = data[vUnit] |
---|
116 |
except: |
---|
117 |
print "Error,",vUnit,"dimension does not exist in",url |
---|
118 |
sys.exit() |
---|
119 |
|
---|
120 |
print |
---|
121 |
print "Subsample factor:",subSample |
---|
122 |
print "Time dim. present:",timeDimPresent |
---|
123 |
|
---|
124 |
# |
---|
125 |
# Check the validity of the X, Y, & T dimensions in the U attrib. vector |
---|
126 |
print "U attibute type:",uUnit |
---|
127 |
a = 0 |
---|
128 |
while 1: |
---|
129 |
try: |
---|
130 |
uVar.dimensions[a] |
---|
131 |
except: |
---|
132 |
break |
---|
133 |
print uUnit,"dimension",a,":",uVar.dimensions[a] |
---|
134 |
if uVar.dimensions[a] not in [xUnit, yUnit, tUnit]: |
---|
135 |
print "Error - " + uVar.dimensions[a] + " is an invalid dimension." |
---|
136 |
sys.exit() |
---|
137 |
a += 1 |
---|
138 |
|
---|
139 |
# Check the validity of the X, Y, & T dimensions in the V attrib. vector |
---|
140 |
print "V attibute type:",vUnit |
---|
141 |
b = 0 |
---|
142 |
while 1: |
---|
143 |
try: |
---|
144 |
vVar.dimensions[b] |
---|
145 |
except: |
---|
146 |
break |
---|
147 |
print vUnit,"dimension",b,":",vVar.dimensions[b] |
---|
148 |
if vVar.dimensions[b] not in [xUnit, yUnit, tUnit]: |
---|
149 |
print "Error - " + vVar.dimensions[b] + " is an invalid dimension." |
---|
150 |
sys.exit() |
---|
151 |
b += 1 |
---|
152 |
|
---|
153 |
# |
---|
154 |
# Establish time step value and time array index location |
---|
155 |
if timeDimPresent: # Deal w/ time dimension if present |
---|
156 |
tArray = tVar[0:tShape] |
---|
157 |
if (tValue <= arcNoData): # No time step specified by user |
---|
158 |
tStep = tShape - 1 # Select single element (last) in time step array |
---|
159 |
tValue = tArray[tShape - 1] |
---|
160 |
else: # Validate user timestep value & set array location |
---|
161 |
tStep = 0 |
---|
162 |
for tElement in tArray: |
---|
163 |
if tElement == tValue: |
---|
164 |
break |
---|
165 |
tStep += 1 |
---|
166 |
else: |
---|
167 |
print "Error: ", tValue, " is an invalid timestep value." |
---|
168 |
print "Valid timestep values: ", tArray |
---|
169 |
sys.exit() |
---|
170 |
print "X array shape = ", xShape |
---|
171 |
print "Y array shape = ", yShape |
---|
172 |
print "time array step = ", tStep |
---|
173 |
print "time step value = ", tValue |
---|
174 |
|
---|
175 |
# |
---|
176 |
# Use time(GMT) in shapefile file name |
---|
177 |
time_handle = strftime('%Y_%m_%d_%H_00_00',(gmtime(tValue))) |
---|
178 |
formatTime = strftime('%m-%d-%Y %H:00 UTC',(gmtime(tValue))) |
---|
179 |
ShapeOut = workspace + "ouba.shp" |
---|
180 |
|
---|
181 |
# |
---|
182 |
# Download X, Y, U, & V data arrays from OPeNDAP server |
---|
183 |
print |
---|
184 |
print "Downloading " + xUnit + ", " + yUnit + ", " + uUnit + ", " + vUnit + " data arrays..." |
---|
185 |
xArray = xVar[0:xShape:subSample] |
---|
186 |
yArray = yVar[0:yShape:subSample] |
---|
187 |
|
---|
188 |
if timeDimPresent: # 3 dimensional vector (T,X,Y) for U and V arrays |
---|
189 |
uArray = uVar[tStep, 0:xShape:subSample, 0:yShape:subSample] |
---|
190 |
vArray = vVar[tStep, 0:xShape:subSample, 0:yShape:subSample] |
---|
191 |
|
---|
192 |
else: # 2 dimensional vector (X,Y) for U and V arrays |
---|
193 |
uArray = uVar[0:xShape:subSample, 0:yShape:subSample] |
---|
194 |
vArray = vVar[0:xShape:subSample, 0:yShape:subSample] |
---|
195 |
|
---|
196 |
print |
---|
197 |
print "Number of Long. values:",xShape |
---|
198 |
print "Number of Lat. values:",yShape |
---|
199 |
print "ASCII temp file:",asciiOut |
---|
200 |
print "Total number of dimensions in",uUnit,"=",a |
---|
201 |
print "Total number of dimensions in",vUnit,"=",b |
---|
202 |
|
---|
203 |
# |
---|
204 |
# Setup temp. process file(s) |
---|
205 |
try: |
---|
206 |
asciiExport = file(asciiOut,'w') |
---|
207 |
except IOError: |
---|
208 |
print "Unable to create temp ASCII data file. Check dir. rights..." |
---|
209 |
sys.exit() |
---|
210 |
|
---|
211 |
# |
---|
212 |
# Load data into ASCII data file (comma delimited -> shapefile - points) |
---|
213 |
x = 0 |
---|
214 |
y = 0 |
---|
215 |
t = 0 |
---|
216 |
# Write column header |
---|
217 |
asciiExport.write("inst,time," + xUnit + "," + yUnit + ",U,V,speed_ms,speed_mph,speed_knots," + "dir,dir_ms\n") |
---|
218 |
# Write X, Y, Z data to temp. ASCII file |
---|
219 |
print |
---|
220 |
print "Writing institution, time, " + xUnit + ", " + yUnit + ", " + uUnit + ", " + vUnit + ", speed(3), and direction(2) data to " + asciiOut + "..." |
---|
221 |
while y < len(yArray): |
---|
222 |
while x < len(xArray): |
---|
223 |
if timeDimPresent: # 3 dimensional vector (T,X,Y) for U and V array |
---|
224 |
asciiExport.write(outPrefix + "," + time_handle + ",") |
---|
225 |
if float(xArray[x]) > 180.: # shift X values from 0 - 360 to -180 - 180 |
---|
226 |
asciiExport.write(str(xArray[x] - 360.) + ",") |
---|
227 |
else: |
---|
228 |
asciiExport.write(str(xArray[x]) + ",") |
---|
229 |
asciiExport.write(str(yArray[y]) + ",") |
---|
230 |
asciiExport.write(str(uArray[t,x,y]) + ",") |
---|
231 |
asciiExport.write(str(vArray[t,x,y]) + ",") |
---|
232 |
# |
---|
233 |
# code to derive Speed from U and V components |
---|
234 |
# (U^2 + V^2)^0.5 |
---|
235 |
speed = (uArray[t,x,y]**2 + vArray[t,x,y]**2)**.5 |
---|
236 |
speed_ms = speed/100 |
---|
237 |
speed_mph = (speed/100) * 2.237 |
---|
238 |
speed_knots = (speed/100) * 1.944 |
---|
239 |
asciiExport.write(str(speed_ms) + "," + str(speed_mph) + "," + str(speed_knots) + ",") |
---|
240 |
# |
---|
241 |
# code to derive Direction from U and V components |
---|
242 |
# atan2(-U,-V) * (180/pi) - 180 |
---|
243 |
direction = atan2(-uArray[t,x,y],-vArray[t,x,y]) * (180/pi) - 180 |
---|
244 |
if direction < 0: |
---|
245 |
azimuth = direction + 360 |
---|
246 |
elif direction > 360: |
---|
247 |
azimuth = direction - 360 |
---|
248 |
ms_azimuth = -azimuth |
---|
249 |
asciiExport.write(str(azimuth) + "," + str(ms_azimuth)) |
---|
250 |
asciiExport.write('\n') |
---|
251 |
|
---|
252 |
else: # 2 dimensional vector (X,Y) for U and V array |
---|
253 |
asciiExport.write(outPrefix + "," + time_handle + ",") |
---|
254 |
if float(xArray[x]) > 180.: # shift X values from 0 - 360 to -180 - 180 |
---|
255 |
asciiExport.write(str(xArray[x] - 360.) + ",") |
---|
256 |
else: |
---|
257 |
asciiExport.write(str(float(xArray[x])) + ",") |
---|
258 |
asciiExport.write(str(yArray[y]) + ",") |
---|
259 |
asciiExport.write(str(uArray[x,y]) + ",") |
---|
260 |
asciiExport.write(str(vArray[x,y])) |
---|
261 |
asciiExport.write('\n') |
---|
262 |
x += 1 |
---|
263 |
x = 0 |
---|
264 |
y += 1 |
---|
265 |
|
---|
266 |
asciiExport.write(outPrefix + "," + formatTime + ",5,10,777,777,NaN,NaN,NaN,NaN,NaN") # bogus data record for TS labeling in MapServer |
---|
267 |
asciiExport.close() |
---|
268 |
|
---|
269 |
# |
---|
270 |
# OGR2OGR Process: ASCII CSV to Shapefile |
---|
271 |
print "Using OGR to create shapefile..." |
---|
272 |
cmd = '/usr/local/bin/ogr2ogr -f "ESRI Shapefile" -where \'U != "' + noData + '"\' ' + ShapeOut + ' /home/jcleary/data/source_files/totals_pydap.vrt' |
---|
273 |
print cmd |
---|
274 |
os.system(cmd) |
---|
275 |
|
---|
276 |
# |
---|
277 |
# MV to NCCOOS folder |
---|
278 |
os.system("mv /home/jcleary/data/codar/seacoos/ouba.shp /home/jcleary/data/codar/nccoos/ouba.shp") |
---|
279 |
os.system("mv /home/jcleary/data/codar/seacoos/ouba.dbf /home/jcleary/data/codar/nccoos/ouba.dbf") |
---|
280 |
os.system("mv /home/jcleary/data/codar/seacoos/ouba.shx /home/jcleary/data/codar/nccoos/ouba.shx") |
---|
281 |
os.system("mv /home/jcleary/data/codar/seacoos/ouba.prj /home/jcleary/data/codar/nccoos/ouba.prj") |
---|
282 |
print |
---|
283 |
print "Move UNC file to be used on NCCOOS.ORG" |
---|
284 |
|
---|
285 |
|
---|
286 |
cleanWorkspace() |
---|
287 |
|
---|
288 |
print |
---|
289 |
print "***" |
---|
290 |
print "*** Data import and conversion complete." |
---|
291 |
print "***" |
---|
292 |
|
---|
293 |
sys.exit() |
---|