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']

"""
Explanation: Plotting with CartoPy and GeoPandas
Converting between GeoPandas and CartoPy for visualizing data.
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.
End of explanation
"""


df.plot()

"""
Explanation: First we'll visualize the map using GeoPandas
End of explanation
"""


# 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()

"""
Explanation: 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.
End of explanation
"""


fig, ax = plt.subplots(subplot_kw={'projection': crs})
ax.add_geometries(df_ae['geometry'], crs=crs)

"""
Explanation: Now that our data is in a CRS based off of CartoPy, we can easily
plot it.
End of explanation
"""


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')

"""
Explanation: Note that we could have easily done this with an EPSG code like so:
End of explanation
"""


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)

"""
Explanation: CartoPy to GeoPandas
Next we'll perform a CRS projection in CartoPy, and then convert it
back into a GeoPandas object.
End of explanation
"""


df_aea = geopandas.GeoDataFrame(df['gdp_pp'], geometry=new_geometries,
                                crs=crs_new.proj4_init)
df_aea.plot()

"""
Explanation: Now that we've created new Shapely objects with the CartoPy CRS,
we can use this to create a GeoDataFrame.
End of explanation
"""


# 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
# Need to provide "zorder" to ensure the points are plotted above the polygons
df_aea_centroids.plot(ax=ax, markersize=5, color='r', zorder=10)

plt.show()

"""
Explanation: 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.
End of explanation
"""


