import numpy as np
import matplotlib.pyplot as plt
import contextily as ctx
import os
import os,sys,inspect
import geopandas as gpd
# I'll start adding two shapefiles for an example. I'll plot the technology district in Buenos Aires polygon and the firms located in the district.
currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) #for relative path reference
distritos_shapes = gpd.read_file(currentdir+'\\data\\distritos-economicos\\distritos_economicos.shp')
distrito_tec=distritos_shapes.loc[distritos_shapes['NOMBRE']=='DISTRITO TECNOLOGICO']
print(distrito_tec.crs)
empresas= gpd.read_file(currentdir+'\\data\\empresas-del-distrito-tecnologico\\empresas-distrito-tecnologico.shp')
print(empresas.crs)
{'proj': 'tmerc', 'lat_0': -34.6297166, 'lon_0': -58.4627, 'k': 0.999998, 'x_0': 100000, 'y_0': 100000, 'ellps': 'intl', 'units': 'm', 'no_defs': True} {'proj': 'tmerc', 'lat_0': -34.6297166, 'lon_0': -58.4627, 'k': 0.999998, 'x_0': 100000, 'y_0': 100000, 'ellps': 'intl', 'units': 'm', 'no_defs': True}
plt.rcParams['figure.figsize'] = [10, 10] #this sets the size of the figure
#Start by defining with matplotlib.plt a new figure and provide a reference name for the axis (in this case myax). Axis are the reference when combining layers
fig, myax = plt.subplots()
#Second, using geopandas' plot function to plot the map in the geopandas dataframe. As passing the plot function,
# add the reference to the named axis. In this case myax
distrito_tec.plot(ax=myax, color='red', alpha=0.2)
#You can add other layers in the same fashion.
empresas.plot(ax=myax, color='blue')
#As new layers are added, the bounds of the axis might change. In order to track the values. The axis() function retrieves the minimum and maximum coordinates of each axis.
myax.axis()
(104170.18040586606, 106817.0328879423, 96749.42815856927, 99937.37243957468)
Now adding the basemap. As explained above, the basemap is collected from the internet. There is an interesting option in contextily to pass a name of a place, similar to a google map seach (see more below), but here we will use the axis coordinates we are already working with in the figure.
#Since contextily only allows passing bounding boxes in both WGS84 (EPSG:4326) and Spheric Mercator (EPSG:3857), it will be necessary to reproject our layers to the appropiate CRS layers first. This is done with the .to_crs method.
distrito_tec_projected= distrito_tec.to_crs(epsg=3857)
empresas_projected= empresas.to_crs(epsg=3857)
# Now repeating the procedure in order to redefine graph with the appropiate coordinates
fig, myax = plt.subplots()
distrito_tec_projected.plot(ax=myax, color='red', alpha=0.2)
empresas_projected.plot(ax=myax, color='blue')
#We will use the add_basemap function featured in geopandas webpage. Notice that the first line retrieves the axis as shown above.
def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
xmin, xmax, ymin, ymax = ax.axis()
basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url)
ax.imshow(basemap, extent=extent, interpolation='bilinear')
# restore original x/y limits
ax.axis((xmin, xmax, ymin, ymax))
add_basemap(myax, zoom=14)
# adding legend with mpatches
import matplotlib.patches as mpatches
red_patch = mpatches.Patch(color='red', label='Distrito Tecnologico')
points_marks = plt.scatter([],[], marker='o',color='blue', label='empresas')
plt.legend(handles=[red_patch,points_marks])
<matplotlib.legend.Legend at 0x1de6192f6a0>