/
nyc_boros.py
72 lines (58 loc) · 2.33 KB
/
nyc_boros.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
"""
Visualizing NYC Boroughs
------------------------
Visualize the Boroughs of New York City with Geopandas.
This example generates many images that are used in the documentation. See
the `Geometric Manipulations <geometric_manipulations>` example for more
details.
First we'll import a dataset containing each borough in New York City. We'll
use the ``datasets`` module to handle this quickly.
"""
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point
from geopandas import GeoSeries, GeoDataFrame
import geopandas as gpd
import geodatasets
np.random.seed(1)
DPI = 100
path_nybb = geodatasets.get_path("nybb")
boros = GeoDataFrame.from_file(path_nybb)
boros = boros.set_index("BoroCode")
boros
##############################################################################
# Next, we'll plot the raw data
ax = boros.plot()
plt.xticks(rotation=90)
plt.savefig("nyc.png", dpi=DPI, bbox_inches="tight")
##############################################################################
# We can easily retrieve the convex hull of each shape. This corresponds to
# the outer edge of the shapes.
boros.geometry.convex_hull.plot()
plt.xticks(rotation=90)
# Grab the limits which we'll use later
xmin, xmax = plt.gca().get_xlim()
ymin, ymax = plt.gca().get_ylim()
plt.savefig("nyc_hull.png", dpi=DPI, bbox_inches="tight")
##############################################################################
# We'll generate some random dots scattered throughout our data, and will
# use them to perform some set operations with our boroughs. We can use
# GeoPandas to perform unions, intersections, etc.
N = 2000 # number of random points
R = 2000 # radius of buffer in feet
# xmin, xmax, ymin, ymax = 900000, 1080000, 120000, 280000
xc = (xmax - xmin) * np.random.random(N) + xmin
yc = (ymax - ymin) * np.random.random(N) + ymin
pts = GeoSeries([Point(x, y) for x, y in zip(xc, yc)])
mp = pts.buffer(R).union_all()
boros_with_holes = boros.geometry - mp
boros_with_holes.plot()
plt.xticks(rotation=90)
plt.savefig("boros_with_holes.png", dpi=DPI, bbox_inches="tight")
##############################################################################
# Finally, we'll show the holes that were taken out of our boroughs.
holes = boros.geometry & mp
holes.plot()
plt.xticks(rotation=90)
plt.savefig("holes.png", dpi=DPI, bbox_inches="tight")
plt.show()