Plotting a map using geopandas and matplotlib

Question

I have small csv that has 6 coordinates from Birmingham England. I read the csv with pandas then transformed it into GeoPandas DataFrame changing my latitude and longitude columns with Shapely Points. I am now trying to plot my GeoDataframe and all I can see are the points. How do I get the Birmingham map represented as well? A good documentation source on GeoPandas would be strongly appreciated too.

from shapely.geometry import Point
import geopandas as gpd
import pandas as pd

df = pd.read_csv('SiteLocation.csv')
df['Coordinates'] = list(zip(df.LONG, df.LAT))
df['Coordinates'] = df['Coordinates'].apply(Point)
# Building the GeoDataframe 
geo_df = gpd.GeoDataFrame(df, geometry='Coordinates')
geo_df.plot()

@Kanak can you add that as an answer?

joris
– joris

2019-01-08 21:50:44 +00:00
Commented Jan 8, 2019 at 21:50 — joris
– joris, Commented Jan 8, 2019 at 21:50
@joris Added as an answer.

keepAlive
– keepAlive

2019-01-10 20:23:31 +00:00
Commented Jan 10, 2019 at 20:23 — keepAlive
– keepAlive, Commented Jan 10, 2019 at 20:23

keepAlive · Accepted Answer · 2019-01-11 15:48:54Z

11

The GeoPandas documentation contains an example on how to add a background to a map (https://geopandas.readthedocs.io/en/latest/gallery/plotting_basemap_background.html), which is explained in more detail below.

You will have to deal with tiles, that are (png) images served through a web server, with a URL like

http://.../Z/X/Y.png, where Z is the zoom level, and X and Y identify the tile

And geopandas's doc shows how to set tiles as backgrounds for your plots, fetching the correct ones and doing all the otherwise difficult job of spatial syncing, etc...

Installation

Assuming GeoPandas is already installed, you need the contextily package in addition. If you are under windows, you may want to pick a look at How to install Contextily?

Use case

Create a python script and define the contextily helper function

import contextily as ctx

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

and play

import matplotlib.pyplot as plt
from shapely.geometry import Point
import geopandas as gpd
import pandas as pd

# Let's define our raw data, whose epsg is 4326
df = pd.DataFrame({
    'LAT'  :[-22.266415, -20.684157],
    'LONG' :[166.452764, 164.956089],
})
df['coords'] = list(zip(df.LONG, df.LAT))

# ... turn them into geodataframe, and convert our
# epsg into 3857, since web map tiles are typically
# provided as such.
geo_df = gpd.GeoDataFrame(
    df, crs  ={'init': 'epsg:4326'},
    geometry = df['coords'].apply(Point)
).to_crs(epsg=3857)

# ... and make the plot
ax = geo_df.plot(
    figsize= (5, 5),
    alpha  = 1
)
add_basemap(ax, zoom=10)
ax.set_axis_off()
plt.title('Kaledonia : From Hienghène to Nouméa')
plt.show()

Note: you can play with the zoom to find the good resolution for the map. E.g./I.e. :

... and such resolutions implicitly call for changing the x/y limits.

edited Jan 11, 2019 at 15:48

answered Jan 8, 2019 at 22:37

keepAlive

6,7055 gold badges29 silver badges43 bronze badges

Sign up to request clarification or add additional context in comments.

8 Comments

keepAlive Over a year ago

Any question @Herc01 ?

joris Over a year ago

Thanks for the answer. What do you mean with the last remark about the zoom level? Normally you should not need to edit the x/y limits manually if changing the zoom level.

keepAlive Over a year ago

@joris, what I meant is that using the value set for the parameter zoom is somehow equivalent to changing the (image) resolution of the map. Do you agree with that?

joris Over a year ago

Yes, that is exactly what the "zoom" means. It might be a bit confusing name as the map itself is not zoomed in (looking at smaller area), but a lower/higher zoom zooms in for the full area and thus resulting in a lower/higher resolution.

joris Over a year ago

sorry, I was planning to add a comment in the edit why I removed it, but forgot that: I know it can certainly be hard to install, but I think here it distracts from the actual answer (and it was also not the only possible way to install, eg for someone using conda that would have been the wrong way). Would you like to make a new question "How to install contextily?", and add your answer on that part there. And then we can link to that question from here (and keep this answer focused on the how).

|

Bruno Vermeulen · Accepted Answer · 2019-11-30 11:04:06Z

Just want to add the use case concerning zooming whereby the basemap is updated according to the new xlim and ylim coordinates. A solution I have come up with is:

First set callbacks on the ax that can detect xlim_changed and ylim_changed
Once both have been detected as changed get the new plot_area calling ax.get_xlim() and ax.get_ylim()
Then clear the ax and re-plot the basemap and any other data

Example for a world map showing the capitals. You notice when you zoom in the resolution of the map is being updated.

import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx


figsize = (12, 10)
osm_url = 'http://tile.stamen.com/terrain/{z}/{x}/{y}.png'
EPSG_OSM = 3857
EPSG_WGS84 = 4326

class MapTools:
    def __init__(self):
        self.cities = gpd.read_file(
            gpd.datasets.get_path('naturalearth_cities'))
        self.cities.crs = EPSG_WGS84
        self.cities = self.convert_to_osm(self.cities)

        self.fig, self.ax = plt.subplots(nrows=1, ncols=1, figsize=figsize)
        self.callbacks_connect()

        # get extent of the map for all cities
        self.cities.plot(ax=self.ax)
        self.plot_area = self.ax.axis()

    def convert_to_osm(self, df):
        return df.to_crs(epsg=EPSG_OSM)

    def callbacks_connect(self):
        self.zoomcallx = self.ax.callbacks.connect(
            'xlim_changed', self.on_limx_change)
        self.zoomcally = self.ax.callbacks.connect(
            'ylim_changed', self.on_limy_change)

        self.x_called = False
        self.y_called = False

    def callbacks_disconnect(self):
        self.ax.callbacks.disconnect(self.zoomcallx)
        self.ax.callbacks.disconnect(self.zoomcally)

    def on_limx_change(self, _):
        self.x_called = True
        if self.y_called:
            self.on_lim_change()

    def on_limy_change(self, _):
        self.y_called = True
        if self.x_called:
            self.on_lim_change()

    def on_lim_change(self):
        xlim = self.ax.get_xlim()
        ylim = self.ax.get_ylim()
        self.plot_area = (*xlim, *ylim)
        self.blit_map()

    def add_base_map_osm(self):
        if abs(self.plot_area[1] - self.plot_area[0]) < 100:
            zoom = 13

        else:
            zoom = 'auto'

        try:
            basemap, extent = ctx.bounds2img(
                self.plot_area[0], self.plot_area[2],
                self.plot_area[1], self.plot_area[3],
                zoom=zoom,
                url=osm_url,)
            self.ax.imshow(basemap, extent=extent, interpolation='bilinear')

        except Exception as e:
            print(f'unable to load map: {e}')

    def blit_map(self):
        self.ax.cla()
        self.callbacks_disconnect()
        cities = self.cities.cx[
            self.plot_area[0]:self.plot_area[1],
            self.plot_area[2]:self.plot_area[3]]
        cities.plot(ax=self.ax, color='red', markersize=3)

        print('*'*80)
        print(self.plot_area)
        print(f'{len(cities)} cities in plot area')

        self.add_base_map_osm()
        self.callbacks_connect()

    @staticmethod
    def show():
        plt.show()


def main():
    map_tools = MapTools()
    map_tools.show()

if __name__ == '__main__':
    main()

Runs on Linux Python3.8 with following pip installs

affine==2.3.0
attrs==19.3.0
autopep8==1.4.4
Cartopy==0.17.0
certifi==2019.11.28
chardet==3.0.4
Click==7.0
click-plugins==1.1.1
cligj==0.5.0
contextily==1.0rc2
cycler==0.10.0
descartes==1.1.0
Fiona==1.8.11
geographiclib==1.50
geopandas==0.6.2
geopy==1.20.0
idna==2.8
joblib==0.14.0
kiwisolver==1.1.0
matplotlib==3.1.2
mercantile==1.1.2
more-itertools==8.0.0
munch==2.5.0
numpy==1.17.4
packaging==19.2
pandas==0.25.3
Pillow==6.2.1
pluggy==0.13.1
py==1.8.0
pycodestyle==2.5.0
pyparsing==2.4.5
pyproj==2.4.1
pyshp==2.1.0
pytest==5.3.1
python-dateutil==2.8.1
pytz==2019.3
rasterio==1.1.1
requests==2.22.0
Rtree==0.9.1
Shapely==1.6.4.post2
six==1.13.0
snuggs==1.4.7
urllib3==1.25.7
wcwidth==0.1.7

Note especially requirement for contextily==1.0rc2

On windows I use Conda (P3.7.3) and don't forget to set the User variables:

GDAL c:\Users\<username>\Anaconda3\envs\<your environment>\Library\share\gdal

PROJLIB c:\Users\<username>\Anaconda3\envs\<your environment>\Library\share

felytic · Accepted Answer · 2019-01-08 13:40:53Z

0

Try df.unary_union. The function will aggregate points into a single geometry. Jupyter Notebook can plot it

answered Jan 8, 2019 at 13:40

felytic

611 silver badge6 bronze badges

1 Comment

Herc01 Over a year ago

Sorry mate -- that was not what I was looking for. But thanks . You look at some other open questions

Collectives™ on Stack Overflow

Plotting a map using geopandas and matplotlib

3 Answers 3

8 Comments

Comments

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

8 Comments

Comments

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related