-
Notifications
You must be signed in to change notification settings - Fork 6
Expand file tree
/
Copy pathpythonvector.py
More file actions
136 lines (99 loc) · 3.97 KB
/
Copy pathpythonvector.py
File metadata and controls
136 lines (99 loc) · 3.97 KB
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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 4 15:16:40 2024
@author: osboxes
"""
from owslib.wfs import WebFeatureService
import geopandas as gpd
import matplotlib.pyplot as plt
'''
# Put the WFS url in a variable
wfsUrl = 'https://service.pdok.nl/cbs/postcode6/2022/wfs/v1_0?request=GetCapabilities&service=WFS'
# Create a WFS object
wfs = WebFeatureService(url=wfsUrl, version='2.0.0')
# Get the title from the object
print(wfs.identification.title)
# Check the contents of the WFS
print(list(wfs.contents))
# Define center point and create bbox for study area
x, y = (173994.1578792833, 444133.60329471016)
xmin, xmax, ymin, ymax = x - 1000, x + 350, y - 1000, y + 350
# Get the features for the study area (using the wfs from the previous code block)
response = wfs.getfeature(typename=list(wfs.contents)[0], bbox=(xmin, ymin, xmax, ymax))
# Save them to disk
with open('data/postal_codes.gml', 'wb') as file:
file.write(response.read())
# Read in again with GeoPandas
pc_gdf = gpd.read_file('data/postal_codes.gml')
# Inspect and plot to get a quick view
print(type(pc_gdf))
pc_gdf.plot()
plt.show()
import json
# Get the WFS of the BAG
wfsUrl = 'https://service.pdok.nl/lv/bag/wfs/v2_0'
wfs = WebFeatureService(url=wfsUrl, version='2.0.0')
layer = list(wfs.contents)[0]
# Define center point and create bbox for study area
x, y = (173994.1578792833, 444133.60329471016)
xmin, xmax, ymin, ymax = x - 500, x + 500, y - 500, y + 500
# Get the features for the study area
# notice that we now get them as json, in contrast to before
response = wfs.getfeature(typename=layer, bbox=(xmin, ymin, xmax, ymax), outputFormat='json')
data = json.loads(response.read())
# Create GeoDataFrame, without saving first
buildingsGDF = gpd.GeoDataFrame.from_features(data['features'])
# Set crs to RD New
buildingsGDF.crs = 28992
# Plot roads and buildings together
pc_layer = pc_gdf.plot(color='grey')
buildingsGDF.plot(ax=pc_layer, color='red')
# Set the limits of the x and y axis
pc_layer.set_xlim(xmin, xmax)
pc_layer.set_ylim(ymin, ymax)
# Save the figure to disk
plt.savefig('./output/postalcoades_roads.png')
# Pandas function that returns the column labels of the DataFrame
print(buildingsGDF.columns)
# Pandas function that returns the first n rows, default n = 5
print(buildingsGDF.head())
# shape area (in the units of the projection)
print(buildingsGDF.area)
# Inspect first
print(buildingsGDF.area > 1000)
# Make the selection, select all rows with area > 1000 m2, and all columns
# Using 'label based' indexing with loc, here with a Boolean array
largeBuildingsGDF = buildingsGDF.loc[buildingsGDF.area > 1000, :]
# Plot
largeBuildingsGDF.plot()
'''
import json
# Define center point and create bbox for study area
x, y = (173994.1578792833, 444133.60329471016)
xmin, xmax, ymin, ymax = x - 500, x + 500, y - 500, y + 500
# Get the WFS of the BAG
wfsUrl = 'https://service.pdok.nl/lv/bag/wfs/v2_0'
wfs = WebFeatureService(url=wfsUrl, version='2.0.0')
layer = list(wfs.contents)[0]
# Get the features for the study area
# notice that we now get them as json, in contrast to before
response = wfs.getfeature(typename=layer, bbox=(xmin, ymin, xmax, ymax), outputFormat='json')
data = json.loads(response.read())
# Create GeoDataFrame, without saving first
buildingsGDF = gpd.GeoDataFrame.from_features(data['features'])
# Set crs to RD New
buildingsGDF.crs = 28992
# Inspect first
print( buildingsGDF['status'] != 'Pand in gebruik' )
# Make the selection, the list of required values can contain more than one item
newBuildingsGDF = buildingsGDF[buildingsGDF['status'] != 'Pand in gebruik']
# Plot the new buildings with a basemap for reference
# based on https://geopandas.org/gallery/plotting_basemap_background.html
import contextily as ctx
# Re-project
newBuildingsGDF = newBuildingsGDF.to_crs(epsg=3857)
# Plot with 50% transparency
ax = newBuildingsGDF.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, zoom=17)
ax.set_axis_off()