/
cartopy_convert.py
105 lines (83 loc) · 3.79 KB
/
cartopy_convert.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
"""
Plotting with CartoPy and GeoPandas
-----------------------------------
Converting between GeoPandas and CartoPy for visualizing data.
`CartoPy <http://scitools.org.uk/cartopy/>`_ is a Python library
that specializes in creating geospatial
visualizations. It has a slightly different way of representing
Coordinate Reference Systems (CRS) as well as constructing plots.
This example steps through a round-trip transfer of data
between GeoPandas and CartoPy.
First we'll load in the data using GeoPandas.
"""
# sphinx_gallery_thumbnail_number = 7
import matplotlib.pyplot as plt
import geopandas
from cartopy import crs as ccrs
path = geopandas.datasets.get_path('naturalearth_lowres')
df = geopandas.read_file(path)
# Add a column we'll use later
df['gdp_pp'] = df['gdp_md_est'] / df['pop_est']
####################################################################
# First we'll visualize the map using GeoPandas
df.plot()
###############################################################################
# Plotting with CartoPy
# =====================
#
# Cartopy also handles Shapely objects well, but it uses a different system for
# CRS. To plot this data with CartoPy, we'll first need to project it into a
# new CRS. We'll use a CRS defined within CartoPy and use the GeoPandas
# ``to_crs`` method to make the transformation.
# Define the CartoPy CRS object.
crs = ccrs.AzimuthalEquidistant()
# This can be converted into a `proj4` string/dict compatible with GeoPandas
crs_proj4 = crs.proj4_init
df_ae = df.to_crs(crs_proj4)
# Here's what the plot looks like in GeoPandas
df_ae.plot()
###############################################################################
# Now that our data is in a CRS based off of CartoPy, we can easily
# plot it.
fig, ax = plt.subplots(subplot_kw={'projection': crs})
ax.add_geometries(df_ae['geometry'], crs=crs)
###############################################################################
# Note that we could have easily done this with an EPSG code like so:
crs_epsg = ccrs.epsg('3857')
df_epsg = df.to_crs(epsg='3857')
# Generate a figure with two axes, one for CartoPy, one for GeoPandas
fig, axs = plt.subplots(1, 2, subplot_kw={'projection': crs_epsg},
figsize=(10, 5))
# Make the CartoPy plot
axs[0].add_geometries(df_epsg['geometry'], crs=crs_epsg,
facecolor='white', edgecolor='black')
# Make the GeoPandas plot
df_epsg.plot(ax=axs[1], color='white', edgecolor='black')
###############################################################################
# CartoPy to GeoPandas
# ====================
#
# Next we'll perform a CRS projection in CartoPy, and then convert it
# back into a GeoPandas object.
crs_new = ccrs.AlbersEqualArea()
new_geometries = [crs_new.project_geometry(ii, src_crs=crs)
for ii in df_ae['geometry'].values]
fig, ax = plt.subplots(subplot_kw={'projection': crs_new})
ax.add_geometries(new_geometries, crs=crs_new)
###############################################################################
# Now that we've created new Shapely objects with the CartoPy CRS,
# we can use this to create a GeoDataFrame.
df_aea = geopandas.GeoDataFrame(df['gdp_pp'], geometry=new_geometries,
crs=crs_new.proj4_init)
df_aea.plot()
###############################################################################
# We can even combine these into the same figure. Here we'll plot the
# shapes of the countries with CartoPy. We'll then calculate the centroid
# of each with GeoPandas and plot it on top.
# Generate a CartoPy figure and add the countries to it
fig, ax = plt.subplots(subplot_kw={'projection': crs_new})
ax.add_geometries(new_geometries, crs=crs_new)
# Calculate centroids and plot
df_aea_centroids = df_aea.geometry.centroid
df_aea_centroids.plot(ax=ax, markersize=5, color='r')
plt.show()