netCDFProcessing/outputDepthStageFromShapeFile.py at master · KevKilo/netCDFProcessing · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#!/usr/bin/env python

#
# avg many shapfiles to make, in this first case, an average May water level
#


import os
import sys

from osgeo import ogr
from osgeo import osr

import numpy
import datetime
import time
import Scientific.IO.NetCDF
from Scientific.IO.NetCDF import NetCDFFile
from scipy import stats

# create list of shapeFiles
#

deBug = False

##shapeFileDir = "/physical/util/shapeFiles/"
shapeFileDir = "/physical/gis/eden/recovery/"

# no. of rows ans columns in the eden grid
#
nrow = 405
ncol = 287

# set for max number of shapefiles to aggregate
# e.g only 31 days in May, but we are looking at 16 Mays from 2000 to 2015

ndays = 500

# create arrrays
#

depth = numpy.zeros( ( 500, 405, 287 ), 'f' )
stage = numpy.zeros( ( 500, 405, 287 ), 'f' )

depthAvg = numpy.zeros( ( 405, 287 ), 'f' )
stageAvg = numpy.zeros( ( 405, 287 ), 'f' )

# get file list
#

##systring = "ls " + shapeFileDir + "eden_epa20000501.shp"
systring = "ls " + shapeFileDir + "*.shp"

print systring

shapeFileList = os.popen( systring, 'r' )

day = 0


for line in shapeFileList.readlines():

        rowList = []
        colList = []

	shapefile = line.rstrip()

        print "Opening %s" % shapefile

	driver = ogr.GetDriverByName("ESRI Shapefile")
	dataSource = driver.Open(shapefile, 0)
	layer = dataSource.GetLayer()

	for feature in layer:
    		currentRow = feature.GetFieldAsInteger( "row" )
    		currentCol = feature.GetFieldAsInteger( "col" )
    		currentWaterDepth = feature.GetFieldAsDouble( "WaterDepth" )
    		currentStage = feature.GetFieldAsDouble( "Stage" )

		rowList.append( currentRow )
		colList.append( currentCol )

		if deBug: print "Row: %d"  % currentRow
		if deBug: print "Cow: %d"  % currentCol
		if deBug: print "Depth: %f"  % currentWaterDepth
		if deBug: print "Stage: %f"  % currentStage

		depth[ day, currentRow, currentCol ] = currentWaterDepth
		stage[ day, currentRow, currentCol ] = currentStage

	day = day + 1

for i in range( 0, day ):
    print "%6.2f" %  depth[ i, 217, 195]