-
Notifications
You must be signed in to change notification settings - Fork 0
/
segments_helper.py
222 lines (187 loc) · 8.01 KB
/
segments_helper.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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 7 12:53:12 2018
@author: willi
"""
from math import sqrt
from datetime import datetime
import numpy as np
from osgeo import ogr, osr
import os
# filedir = r'D:\MASTER_GEOTECH_2S\LastCourse\FinalProject\Data'
# #path to the shapefile
# points_shp_filename = os.path.join(filedir, 'shp','participant_6_pr.shp')
# # Path to output shapefile
# segments_shp_filename = os.path.join(filedir, 'shp','segments_participant_6.shp')
def generate_segments(points_shp_filename, segments_shp_filename):
# Load the shapefile and extract the layer
driver = ogr.GetDriverByName('ESRI Shapefile')
shapefile = driver.Open(points_shp_filename, 0)
if shapefile == None:
print(f"Shapefile '{points_shp_filename}' contains no data. Skipping...")
return
layer = shapefile.GetLayer()
# If the output file already exists, delete it
if os.path.exists(segments_shp_filename):
print(f"A file with the name '{segments_shp_filename}' already exists. Deleting...")
driver.DeleteDataSource(segments_shp_filename)
# Create the output shapefile
lines_ds = driver.CreateDataSource(segments_shp_filename)
# Set the spatial reference, WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(32632)
#=================
# Lines
#=================
# Create a new layer within the output shapefile
lines = lines_ds.CreateLayer('lines', srs, ogr.wkbLineString)
fields = [
{'name': 'id_person', 'type': ogr.OFTString},
{'name': 'start_time', 'type': ogr.OFTString},
{'name': 'end_time', 'type': ogr.OFTString},
{'name': 'distance', 'type': ogr.OFTReal},
{'name': 'duration', 'type': ogr.OFTReal},
{'name': 'speed', 'type': ogr.OFTReal},
{'name': 'accuracy', 'type': ogr.OFTReal},
{'name': 'activity', 'type': ogr.OFTString},
{'name': 'confidence', 'type': ogr.OFTReal},
{'name': 'tmax', 'type': ogr.OFTReal},
{'name': 'tmin', 'type': ogr.OFTReal},
{'name': 'precip', 'type': ogr.OFTReal},
]
for field in fields:
field_definition = ogr.FieldDefn(field['name'], ogr.OFTReal['type'])
lines.CreateField(field_definition)
#===========================================
# Function to calculate distance
#===========================================
def euclideandist(x1,y1,x2,y2):
deltax = x1 - x2
deltay = y1 - y2
d = sqrt(deltax**2 + deltay**2)
return(d)
# Secondly, we wrote a function that calculates the duration of each segment. This is simply the difference between the timestamp of the startpoint of the segment and the timestamp of the endpoint of the segment. In our data, the timestamp of each point is stored as a string. The module datetime was used to transform these strings into datetime objects, such that substractions of timestamps coudl be done correctly. The returned value will be duration in seconds.
# Define the duration function
def timedif(time1, time2):
# Convert to datetime objects
time1 = datetime.strptime(time1, '%Y-%m-%d %H:%M:%S')
time2 = datetime.strptime(time2, '%Y-%m-%d %H:%M:%S')
# Convert to numeric (seconds from 01-01-1970 UTC)
time1 = time1.timestamp()
time2 = time2.timestamp()
# Calculate time difference
dT = time2 - time1
return dT
#======================
# Get id of participants
#======================
idn_list = []
for feat in layer:
idn = feat.GetField("id_person")
idn_list.append(idn)
layer.ResetReading()
idents = np.unique(idn_list)
#======================
# For loop segments
#======================
# Select an individual goose base on its unique ID and create the corresponding segments, including attributes
for i in idents:
# Create empty lists to store the geometry and attribute values of each point in
x_list = [] # X coordinates
y_list = [] # Y coordinates
time_list = [] # Timestamps
accur_list = [] # Accuracy
activ_list = [] # Activity
confi_list = [] # Confidence
Tmax_list = [] # Tmeperature max
Tmin_list = [] # Temperature min
prec_list = [] # Precipitation
# Store the geometry and attribute values of each point in lists
for ft in layer:
idn = ft.GetField('id_person')
if idn == i:
pt = ft.geometry()
x = pt.GetX()
y = pt.GetY()
time= ft.GetField('timestamp')
accur = ft.GetField('accuracy')
activ = ft.GetField('activity')
confi = ft.GetField('confidence')
Tmax = ft.GetField('tmax')
Tmin = ft.GetField('tmin')
prec = ft.GetField('precip')
x_list.append(x)
y_list.append(y)
time_list.append(time)
accur_list.append(accur)
activ_list.append(activ)
confi_list.append(confi)
Tmax_list.append(Tmax)
Tmin_list.append(Tmin)
prec_list.append(prec)
layer.ResetReading()
# Create each segment one by one
for j in range(0, len(x_list)-1):
# Get geometry and attribute values of startpoint and endpoint of the segment
stx = x_list[j]
enx = x_list[j+1]
sty = y_list[j]
eny = y_list[j+1]
start_time = time_list[j]
end_time = time_list[j+1]
stTmax = Tmax_list[j]
enTmax = Tmax_list[j+1]
stTmin = Tmin_list[j]
enTmin = Tmin_list[j+1]
stprec = prec_list[j]
enprec = prec_list[j+1]
# Create a segment between the start and end point
seg = ogr.Geometry(ogr.wkbLineString)
seg.AddPoint(stx, sty)
seg.AddPoint(enx, eny)
# Create a new feature
feat = ogr.Feature(lines.GetLayerDefn())
# Add the linestring geometry to the feature
feat.SetGeometry(seg)
# Add the corresponding goose ID to the 'idn_ident' field of the feature
feat.SetField('id_person', i)
# Add start and end time by segment
feat.SetField('start_time', start_time)
feat.SetField('end_time', end_time)
# Calculate distance and add this value to the corresponding field
distance = euclideandist(stx, sty, enx, eny)
feat.SetField('distance', distance)
# Calculate duration and add this value to the corresponding field
sttime1 = start_time.split('.')[0]
entime1 = end_time.split('.')[0]
duration = timedif(sttime1, entime1)
feat.SetField('duration', duration)
# Activity, I will select the start point
accurl = accur_list[j]
feat.SetField('accuracy', accurl)
# Activity, I will select the start point
act = activ_list[j]
feat.SetField('activity', act)
# Confidence, I will select the start point
conf = confi_list[j]
feat.SetField('confidence', conf)
# Calculate speed and add this value to the corresponding field
if duration == 0:
speed = 0
else:
speed = distance/duration
feat.SetField('speed', speed)
# Calculate average weather values and add these values to the corresponding fields
Tmaxi = (stTmax + enTmax)/2
feat.SetField('tmax', Tmaxi)
Tmini = (stTmin + enTmin)/2
feat.SetField('tmin', Tmini)
prec = (stprec + enprec)/2
feat.SetField('precip', prec)
# Create the feature with its field
lines.CreateFeature(feat)
# Dereference the feature
feat = None
# Close the data source
lines_ds = None
print('Done!')