Exercise outline
We want to create a station class (data structure) that can hold the data of one station. Then we make a list containing all the stations in the NetCdf file and fill it with the data of the NetCdf file.
Create list of stations
Add the following import statement to the region Import libs:
# DeltaShell libs from SharpMap.CoordinateSystems.Transformations import GeometryTransform
The GeometryTransform class is needed to transform coordinates from one coordinatesystem to another coordinatesystem.
And add this code to the end of the script :
#region Create station list
class Station:
# Station definition (name, geometry and componentInformation (amplitudes and phases))
def __init__(self, name, point, amplitudeArray, phaseArray):
self.name = name.strip()
self.geometry = point
self.amplitudeArray = amplitudeArray
self.phaseArray = phaseArray
# Check for coordinate system
stationCoordinateSystem = file.GetAttributeValue(file.GetVariableByName("crs"),"coord_ref_sys_name")
print "Station information is in " + stationCoordinateSystem
print "Model information is in " + fmModel.CoordinateSystem.Name
# Create transformation if needed
if (stationCoordinateSystem != fmModel.CoordinateSystem.Name):
wgs84 = Map.CoordinateSystemFactory.CreateFromEPSG(4326) # coordinateSystem of stations (wgs 84)
transformation = Map.CoordinateSystemFactory.CreateTransformation(wgs84, fmModel.CoordinateSystem).MathTransform
else:
transformation = None
stationList = []
for i in range(nrOfStations):
phaseComponents = GetComponentData(phaseData, i, nrOfComponents)
amplitudeComponents = GetComponentData(amplitudeData, i, nrOfComponents)
geometry = Point(lonData[i],latData[i])
if (transformation != None):
geometry = GeometryTransform.TransformPoint(geometry, transformation)
# create station
station = Station(stationNames[i], geometry, amplitudeComponents, phaseComponents)
stationList.append(station)
print "Number of stations found : " + str(len(stationList))
#endregion
The following code declares a station class (definition of a station object). It states that a station consists of the following properties : name, geometry, amplitudeArray and phaseArray
class Station: # Station definition (name, geometry and componentInformation (amplitudes and phases)) def __init__(self, name, point, amplitudeArray, phaseArray): self.name = name.strip() self.geometry = point self.amplitudeArray = amplitudeArray self.phaseArray = phaseArray
The function _init_ is called when you create an object of type station (constructor). In this case you can only create a station when you specify all the properties.
After declaring the Station class we continue by comparing the coordinatesystem of the NetCdf file and fmModel :
# Check for coordinate system
stationCoordinateSystem = file.GetAttributeValue(file.GetVariableByName("crs"),"coord_ref_sys_name")
print "Station information is in " + stationCoordinateSystem
print "Model information is in " + fmModel.CoordinateSystem.Name
If the coordinatesystem differs (when using harlingen model for instance (Rd new)) we create transformation between the two coordinatesystem :
# Create transformation if needed if (stationCoordinateSystem != fmModel.CoordinateSystem.Name): wgs84 = Map.CoordinateSystemFactory.CreateFromEPSG(4326) # coordinateSystem of stations (wgs 84) transformation = Map.CoordinateSystemFactory.CreateTransformation(wgs84, fmModel.CoordinateSystem).MathTransform else: transformation = None
Now we can create stations and add them to the "stationList" with the following code :
stationList = [] for i in range(nrOfStations): phaseComponents = GetComponentData(phaseData, i, nrOfComponents) amplitudeComponents = GetComponentData(amplitudeData, i, nrOfComponents) geometry = Point(lonData[i],latData[i]) if (transformation != None): geometry = GeometryTransform.TransformPoint(geometry, transformation) # create station station = Station(stationNames[i], geometry, amplitudeComponents, phaseComponents) stationList.append(station) print "Number of stations found : " + str(len(stationList))
