-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmod17_Parser.py
175 lines (162 loc) · 7.54 KB
/
mod17_Parser.py
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
import gdal
import gdalconst
import pandas as pd
import pdb
class BilFile(object):
# This is the class used to open a .bil file and make georeferenced data
# available. Raster values for a coordinate are extracted using the
# extract_coord_val method
def __init__(self, bil_file):
# Construct the BilFile object and initialize its
# properties
self.bil_file = bil_file
self.hdr_file = bil_file.split('.')[0]+'.hdr'
# Register EHdr driver needed for these files - automatic in python?
# gdal.GetDriverByName('EHdr').Register()
# Open the file and put in dataset
dataset = gdal.Open(self.bil_file, gdalconst.GA_ReadOnly)
band = dataset.GetRasterBand(1)
self.ncol = dataset.RasterXSize
self.nrow = dataset.RasterYSize
# Get georeferencing information and assign to object
# properties
geotransform = dataset.GetGeoTransform()
self.originX = geotransform[0]
self.originY = geotransform[3]
self.pixelWidth = geotransform[1]
self.pixelHeight = geotransform[5]
self.data = band.ReadAsArray()
def extract_coord_val(self, lat, lon):
# Method for extracting raster values at given coordinates
y = int((lat - self.originY)/self.pixelHeight)
x = int((lon - self.originX)/self.pixelWidth)
return self.data[y, x]
# Function for extracting daily PRISM data
def getDailyPrism(year, metdata, data_path, coords_file):
# Read in site coordinates, get date range and create a DataFrame
#to fill
pnts = pd.read_csv(coords_file)
drange = pd.date_range('1-1-{0}'.format(year),
'12-31-{0}'.format(year), freq='D')
df = pd.DataFrame(index=drange, columns=pnts.sitecode)
for i in range(len(drange)):
# Create a tuple to fill the file dates in,
# pad month & day with zeros
ymd_tuple = (str(drange.year[i]),
str(drange.month[i]).zfill(2),
str(drange.day[i]).zfill(2))
# Create GDAL vsizip file path (read directly from zip archive)
# See https://trac.osgeo.org/gdal/wiki/UserDocs/ReadInZip
# If for some reason this vsizip interface won't work, extract the
# archive and remove the 'vsizip' part of pathname
bil_file = (r'/vsizip/' + data_path +
r'PRISM_{0}_stable_4kmD2_{1}0101_{1}1231_bil.zip/'.format(
metdata, year) +
r'PRISM_{0}_stable_4kmD2_{1}{2}{3}_bil.bil'.format(metdata, *ymd_tuple))
bil_ds = BilFile(bil_file)
for j in range(len(pnts.index)):
pt_val = bil_ds.extract_coord_val(pnts.lat[j], pnts.lon[j])
df.iloc[i, j] = pt_val
return df
# Function for extracting daily PRISM data from provisional files
def getDailyPrismProvis(year, month, metdata, data_path, bil_name, coords_file):
# Read in site coordinates, get date range and create a DataFrame
#to fill
pnts = pd.read_csv(coords_file)
drange = pd.date_range('{0}-1-1'.format(year),
'{0}-12-31'.format(year), freq='D')
drange = drange[drange.month==month]
df = pd.DataFrame(index=drange, columns=pnts.sitecode)
for i in range(len(drange)):
# Create a tuple to fill the file dates in,
# pad month & day with zeros
ymd_tuple = (str(drange.year[i]),
str(drange.month[i]).zfill(2),
str(drange.day[i]).zfill(2))
# Create GDAL vsizip file path (read directly from zip archive)
# See https://trac.osgeo.org/gdal/wiki/UserDocs/ReadInZip
# If for some reason this vsizip interface won't work, extract the
# archive and remove the 'vsizip' part of pathname
bil_file = (r'/vsizip/' + data_path + bil_name + '/' +
r'PRISM_{0}_provisional_4kmD2_{1}{2}{3}_bil.bil'.format(metdata,
*ymd_tuple))
bil_ds = BilFile(bil_file)
for j in range(len(pnts.index)):
pt_val = bil_ds.extract_coord_val(pnts.lat[j], pnts.lon[j])
df.iloc[i, j] = pt_val
return df
# Function for extracting monthly PRISM data
# Note that this uses 1981-2015 files and will need to be altered if different
# files are used
def getMonthlyPrism( metdata, data_path, coords_file ):
# Read in site coordinates, get date range and create a DataFrame
#to fill
pnts = pd.read_csv(coords_file)
drange = pd.date_range('1-1-1981', '9-30-2015', freq='M')
df = pd.DataFrame(index=drange, columns=pnts.sitecode)
for i in range(len(drange)):
# Create a tuple to fill the file dates in,
# pad month & day with zeros
ym_tuple = (str(drange.year[i]),
str(drange.month[i]).zfill(2))
# Create GDAL vsizip file path (read directly from zip archive)
# See https://trac.osgeo.org/gdal/wiki/UserDocs/ReadInZip
# If for some reason this vsizip interface won't work, extract the
# archive and remove the 'vsizip' part of pathname
if metdata=='ppt':
bil_file = (r'/vsizip/' + data_path +
r'PRISM_{0}_stable_4kmM3_198101_201509_bil.zip/'.format(
metdata) +
r'PRISM_{0}_stable_4kmM3_{1}{2}_bil.bil'.format(
metdata, *ym_tuple))
elif metdata=='tmean':
bil_file = (r'/vsizip/' + data_path +
r'PRISM_{0}_stable_4kmM2_198101_201509_bil.zip/'.format(
metdata) +
r'PRISM_{0}_stable_4kmM2_{1}{2}_bil.bil'.format(
metdata, *ym_tuple))
bil_ds = BilFile(bil_file)
for j in range(len(pnts.index)):
pt_val = bil_ds.extract_coord_val(pnts.lat[j], pnts.lon[j])
df.iloc[i, j] = pt_val
return df
# Function for extracting monthly PRISM data from provisional files
def getMonthlyPrismProvis(year, metdata, data_path, bil_name, coords_file):
# Read in site coordinates, get date range and create a DataFrame
#to fill
pnts = pd.read_csv(coords_file)
drange = pd.date_range('{0}-10-01'.format(year),
'{0}-12-31'.format(year), freq='M')
#drange = drange[drange.month==month]
df = pd.DataFrame(index=drange, columns=pnts.sitecode)
for i in range(len(drange)):
# Create a tuple to fill the file dates in,
# pad month & day with zeros
ym_tuple = (str(drange.year[i]),
str(drange.month[i]).zfill(2))
# Create GDAL vsizip file path (read directly from zip archive)
# See https://trac.osgeo.org/gdal/wiki/UserDocs/ReadInZip
# If for some reason this vsizip interface won't work, extract the
# archive and remove the 'vsizip' part of pathname
bil_file = (r'/vsizip/' + data_path + bil_name + '/' +
r'PRISM_{0}_provisional_4kmM2_{1}{2}_bil.bil'.format(metdata,
*ym_tuple))
#pdb.set_trace()
bil_ds = BilFile(bil_file)
for j in range(len(pnts.index)):
pt_val = bil_ds.extract_coord_val(pnts.lat[j], pnts.lon[j])
df.iloc[i, j] = pt_val
return df
# Function for extracting 30 year normal PRISM precip data
def get30yrPrismPrecip(data_path, coords_file):
# Read in site coordinates, get date range and create a DataFrame
#to fill
pnts = pd.read_csv(coords_file)
df = pd.DataFrame(index=pnts.sitecode, columns=['30yr_ppt'])
bil_file = (data_path +
r'PRISM_ppt_30yr_normal_800mM2_annual_bil.bil')
bil_ds = BilFile(bil_file)
for j in range(len(pnts.index)):
precip = bil_ds.extract_coord_val(pnts.lat[j], pnts.lon[j])
df.loc[pnts.sitecode[j], '30yr_ppt'] = precip
return df