Commit bf193538 authored by Gerrit Erichsen's avatar Gerrit Erichsen

Differentiated between types in evaluation and adjusted modeled data to measured data types

parent ca184d54
This diff is collapsed.
......@@ -17,6 +17,7 @@ class WeatherStation:
m_lon = ""
m_measuredData = []
m_modelData = []
m_minutes = []
m_mae = 0.
m_mbe = 0.
m_rmse = 0.
......@@ -24,12 +25,24 @@ class WeatherStation:
m_obsMin = 0.
m_obsAvg = 0.
m_error = 0.2
m_type = 'temperature'
def __init__(self, name, identifier, lat, lon):
def __init__(self, name, identifier, lat, lon, dataType = 'temperature'):
self.m_name = name
self.m_id = identifier
self.m_lat = lat
self.m_lon = lon
self.m_measuredData = []
self.m_modelData = []
self.m_minutes = []
self.m_mae = 0.
self.m_mbe = 0.
self.m_rmse = 0.
self.m_obsMax = 0.
self.m_obsMin = 0.
self.m_obsAvg = 0.
self.m_error = 0.2
self.m_type = dataType
def getName(self):
return self.m_name
......@@ -50,9 +63,17 @@ class WeatherStation:
def getModelData(self):
data = []
for element in self.m_modelData:
#TODO: make type dependent
data.append(element - 273.15)
if self.m_type == 'temperature':
data.append(element - 273.15)
else:
data.append(element)
return data
def getModelDataLen(self):
return len(self.m_modelData)
def getModelDataAtIndex(self, index):
if index >= 0 and index < len(self.m_modelData):
return self.m_modelData[index]
return 0.
def getDifferences(self):
differences = []
i = 0
......@@ -60,23 +81,111 @@ class WeatherStation:
if (len(self.m_modelData) < totalLength):
totalLength = len(self.m_modelData)
while (i < totalLength):
#TODO: make type dependent
if (self.m_measuredData[i] > -500. and
self.m_measuredData[i] < 2000.):
differences.append(self.m_measuredData[i] - self.m_modelData[i] \
+ 273.15)
value = self.m_measuredData[i] - self.m_modelData[i]
if (self.m_type == 'temperature'):
value += 273.15
elif (self.m_type == 'windDirection'):
if value > 180. :
value -= 180.
elif value < -180. :
value += 180.
elif self.m_type == 'solar' \
or self.m_type == 'diffuseIrridiation' \
or self.m_type == 'globalIrridiation':
joulePerCm2ToWhPerM2 = 1 * 10000 / 3600
observed = self.m_measuredData[i] * joulePerCm2ToWhPerM2
modeled = 0.
if i > 0 :
modeled = self.m_modelData[i] * self.m_minutes[i] / 60\
+ self.m_modelData[i - 1] \
* (60 - self.m_minutes[i - 1]) / 60
value = observed - modeled
differences.append(value)
i += 1
return differences;
def getHeightIndices(height):
""""(float) -> (str, str)
function returns indices of COSMO-DE's main planes, that lie
around the specified height"""
lower = ''
upper = ''
if height <= 10.:
lower = '50'
elif height < 35.72:
lower = '50'
upper = '49'
elif height == 35.72:
lower = '49'
elif height < 73.03:
lower = '49'
upper = '48'
elif height == 73.03:
lower = '48'
elif height < 122.32:
lower = '48'
upper = '47'
elif height == 122.32:
lower = '47'
elif height < 183.93:
lower = '47'
upper = '46'
elif height == 183.93:
lower = '46'
elif height < 258.21:
lower = '46'
upper = '45'
elif height == 258.21:
lower = '45'
elif height < 345.53:
lower = '45'
upper = '44'
elif height >= 345.53:
lower = '44'
return lower, upper
def getHeightFromIndex(index):
""""(str) -> (float)
function returns height of COSMO-DE's main plane with index index"""
if index == '50':
return 10.
elif index == '49':
return 35.72
elif index == '48':
return 73.03
elif index == '47':
return 122.32
elif index == '46':
return 183.93
elif index == '45':
return 258.21
else:
return 345.53
def setModelData(self, data):
if len(data) > 0:
self.m_modelData = data
def setType(self, dataType):
self.m_type = dataType
def readInCdcFile(self, fileName, targetYear, column):
self.m_measureData = [] #clear history -> expercience shows to be necessary
file = open(fileName, "r")
content = file.readlines()
file.close()
lastHour = 23 #entry that was before current one, needed for missing data
counter = 0
for line in content:
entries = line.split(';')
if entries[1].startswith(targetYear):
hour = int(entries[1][-2:])
if self.m_type == 'solar'\
or self.m_type == 'diffuseIrridiation' \
or self.m_type == 'globalIrridiation':
hour = int(entries[1][-5:-3])
self.m_minutes.append(int(entries[1][-2:]))
else:
hour = int(entries[1][-2:])
if hour > lastHour + 1:
i = 0
diff = hour - lastHour - 1
......@@ -91,11 +200,31 @@ class WeatherStation:
self.m_measuredData.append(-999.)
i += 1
self.m_measuredData.append(float(entries[column]))
if self.m_type == 'windDirection' and counter == 0:
print(self.m_type, column, entries[column], self.m_measuredData)
lastHour = hour
counter += 1
def readInModelFile(self, h5FileName, dataType, row, col):
with h5py.File(h5FileName, 'r') as h5File:
self.m_modelData = h5File[dataType][:,row,col]
def readInModelFile(self, h5FileName, dataType, row, col, \
h5FileName2 = '', height = 10.):
if self.m_name == 'Brocken':
print('Brocken, alter', h5FileName, h5FileName2, height)
if h5FileName2 == '':
with h5py.File(h5FileName, 'r') as h5File:
self.m_modelData = h5File[dataType][:,row,col]
else:
data1 = []
data2 = []
height1 = WeatherStation.getHeightFromIndex(h5FileName[-5:-3])
height2 = WeatherStation.getHeightFromIndex(h5FileName2[-5:-3])
with h5py.File(h5FileName, 'r') as h5File1:
data1 = h5File1[dataType][:,row,col]
with h5py.File(h5FileName2, 'r') as h5File2:
data2 = h5File2[dataType][:,row,col]
for lower, upper in zip(data1, data2):
self.m_modelData.append(np.interp(height, [height1, height2], \
[lower, upper]))
def calculateStats(self):
differences = self.getDifferences();
......@@ -106,6 +235,11 @@ class WeatherStation:
for element in self.m_measuredData:
if (element > -500. and element < 2000.):
values.append(element);
if self.m_type == 'solar' \
or self.m_type == 'diffuseIrridiation' \
or self.m_type == 'globalIrridiation':
joulePerCm2ToWhPerM2 = 1 * 10000 / 3600
values[-1] *= joulePerCm2ToWhPerM2
valueArray = np.array(values);
self.m_obsMin = valueArray.min()
self.m_obsMax = valueArray.max()
......@@ -133,8 +267,21 @@ class WeatherStation:
if (self.m_measuredData[i] > -500. \
and self.m_measuredData[i] < 2000.):
valuesObs.append(self.m_measuredData[i])
valuesObsBounds.append(self.m_error)
valuesMod.append(self.m_modelData[i] - 273.15)
valuesObsBounds.append(self.m_error)
valueMod = self.m_modelData[i]
if self.m_type == 'temperature':
valueMod -= 273.15
elif self.m_type == 'solar' \
or self.m_type == 'diffuseIrridiation' \
or self.m_type == 'globalIrridiation':
joulePerCm2ToWhPerM2 = 1 * 10000 / 3600
valuesObs[-1] *= joulePerCm2ToWhPerM2
valueMod = 0.
if i > 0 :
valueMod= self.m_modelData[i] * self.m_minutes[i] / 60\
+ self.m_modelData[i - 1] \
* (60 - self.m_minutes[i - 1]) / 60
valuesMod.append(valueMod)
xAxis.append(i)
i += 1
plt.figure()
......@@ -144,7 +291,7 @@ class WeatherStation:
upper = np.array(valuesObs) + np.array(valuesObsBounds)
lower = np.array(valuesObs) - np.array(valuesObsBounds)
plt.fill_between(xAxis, lower, upper)
plt.title(self.getName())
plt.title(self.getName() + ' ' + self.m_type + ' Run')
plt.legend()
plt.show()
plt.figure()
......@@ -159,4 +306,5 @@ class WeatherStation:
self.m_obsMax - self.m_obsMin, \
self.m_obsAvg)
plt.annotate(annotation, (maxX / 2,0))
plt.title(self.getName() + ' ' + self.m_type + ' Scatter')
plt.show()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment