-
Notifications
You must be signed in to change notification settings - Fork 0
/
shapefile_helper.py
153 lines (134 loc) · 5.47 KB
/
shapefile_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
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 6 16:21:21 2018
@author: willi
"""
import os
import csv
from osgeo import ogr, osr
def json_to_shapefile(data_points, output_filename, identifier=None):
if identifier == None:
head, tail = os.path.split(output_filename)
identifier = tail.replace(".shp", "")
# Initialize OGR driver
driver = ogr.GetDriverByName('Esri Shapefile')
# If the output file already exists, delete it
if os.path.exists(output_filename):
print(f"A file with the name '{output_filename}' already exists. Deleting...")
driver.DeleteDataSource(output_filename)
# Create the shapefile
print("Creating shapefile...")
shapefile = driver.CreateDataSource(output_filename)
# Create main layer with WGS 84 reference system
spatial_reference_system = osr.SpatialReference()
spatial_reference_system.ImportFromEPSG(4326)
layer = shapefile.CreateLayer('locations', spatial_reference_system, ogr.wkbPoint)
# Add fields to the layer
field_definitions = [
ogr.FieldDefn('id', ogr.OFTString),
ogr.FieldDefn('id_person', ogr.OFTString),
ogr.FieldDefn('timestamp', ogr.OFTString),
ogr.FieldDefn('latitude', ogr.OFTReal),
ogr.FieldDefn('longitude', ogr.OFTReal),
ogr.FieldDefn('accuracy', ogr.OFTReal),
# ogr.FieldDefn('velocity', ogr.OFTReal),
# ogr.FieldDefn('heading', ogr.OFTReal),
# ogr.FieldDefn('altitude', ogr.OFTReal),
# ogr.FieldDefn('v_accuracy', ogr.OFTReal),
ogr.FieldDefn('activity', ogr.OFTString),
ogr.FieldDefn('confidence', ogr.OFTReal),
]
for field_definition in field_definitions:
layer.CreateField(field_definition)
# Add a feature to the layer for each data_point
print("Adding features...")
for index, data_point in enumerate(data_points):
# Initialize a new feature (geometry)
feature = ogr.Feature(layer.GetLayerDefn())
# Set the geometry of the new feature
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(float(data_point["lon"]), float(data_point["lat"]))
feature.SetGeometry(point)
# Set the attributes of the new feature
# TODO: Handle "" values for OFTReal fields (should not become zeros)
feature.SetField('id', index)
feature.SetField('id_person', identifier)
feature.SetField('timestamp', str(data_point["timestamp"]))
feature.SetField('latitude', data_point["lat"])
feature.SetField('longitude', data_point["lon"])
feature.SetField('accuracy', data_point["accuracy"])
# feature.SetField('velocity', data_point["velocity"])
# feature.SetField('heading', data_point["heading"])
# feature.SetField('altitude', data_point["altitude"])
# feature.SetField('v_accuracy', data_point["vertical_accuracy"])
feature.SetField('activity', data_point["activity"])
feature.SetField('confidence', data_point["activity_confidence"])
# Add the feature to the layer
layer.CreateFeature(feature)
def csv_to_shp(filename_csv, name_person, filename_out_shp):
#============================
#Importing shapefile
#============================
infile = open(filename_csv, 'r') # CSV file
table = []
csv_reader = csv.reader(infile)
for row in csv_reader:
table.append(row)
infile.close()
#table = table[0:1000]
#len(table)
#============================
#Creating shapefile
#============================
# Now convert it to a shapefile with OGR
driver = ogr.GetDriverByName('Esri Shapefile')
#if the output shapefile that we create already exists, delete it
if os.path.exists(filename_out_shp):
print('deleting')
driver.DeleteDataSource(filename_out_shp)
ds = driver.CreateDataSource(filename_out_shp)
#reference system
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
layer = ds.CreateLayer('pt', srs , ogr.wkbPoint)
# Add new fields
id_def = ogr.FieldDefn('id', ogr.OFTString) #ID
layer.CreateField(id_def)
idp_def = ogr.FieldDefn('id_person', ogr.OFTString) #ID
layer.CreateField(idp_def)
tim_def = ogr.FieldDefn('timestamp', ogr.OFTString) #timestamp
layer.CreateField(tim_def)
lat_def = ogr.FieldDefn('latitud', ogr.OFTReal) #latitud
layer.CreateField(lat_def)
lon_def = ogr.FieldDefn('longitud', ogr.OFTReal) #latitud
layer.CreateField(lon_def)
act_def = ogr.FieldDefn('activity', ogr.OFTString) #latitud
layer.CreateField(act_def)
#for loop geometries
j=0
for rowa in table:
if j==0:
j = j+1
continue
else:
# Create a new feature (attribute and geometry)
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(float(rowa[2]), float(rowa[1]))
feat = ogr.Feature(layer.GetLayerDefn())
feat.SetGeometry(point)
# add ID
feat.SetField('id', j)
# add ID
feat.SetField('id_person', name_person)
# add timestamp
feat.SetField('timestamp', rowa[0])
# add latitud
feat.SetField('latitud', float(rowa[1]))
# add longitud
feat.SetField('longitud', float(rowa[2]))
# add longitud
feat.SetField('activity', rowa[-2])
# Make a geometry, from Shapely object
layer.CreateFeature(feat)
feat = point = None # destroy these
j = j+1