Geospatial data mapping with python

  • Post author:
  • Post category:GIS
  • Post comments:6 Comments
Geospatial data visualisation and mapping with python
Geospatial data visualisation and mapping with python

Introduction

Mapping of geospatial data with python is really fun to work. The main library for this work is GDAL library. Perhaps, directly working with GDAL is bit complicated. we will mostly use geopandas of python. Geopandas mostly act as a wrapper. Under the hood it uses GDAL library. It also uses shapely and other GIS related packages. In addition, geopandas can handle the pandas dataframe as well.

Importing Packages

At first we will import all the packages to make a geographical map. Geopandas is the main package to make map of geospatial data.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Point
import os
import matplotlib.gridspec as gridspec
import matplotlib.patheffects as PathEffects
from scipy.interpolate import griddata
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

Making point shape file from geospatial data

We will make point shape file from coordinate data. There is a parameter in each location. We will change colour of the point shape file based on that value later on. Point shapefile is saved in a variable named ‘sites’.

binpath_b = '/Users/ihasan/Downloads/Blogging/Mapping in python with geopandas/'
fn_b = os.path.join(binpath_b,'point_inf.csv')
df = pd.read_csv(fn_b)
df['Coordinates'] = list(zip(df.Long, df.Lat))
df['Coordinates'] = df['Coordinates'].apply(Point)
sites = gpd.GeoDataFrame(df, geometry='Coordinates')

Loading other shapefile directory

Addressing the geospatial data directory of the district of Bangladesh and the Forest shapefile as well. The district shapefile of Bangladesh is saved in the variable named ‘Bangladesh’ and the forest shapefile is saved in variable named ‘fr’.

sdir1 = os.path.dirname('/Users/ihasan/Downloads/Rivers of Bangladesh/Bangladesh/BGD_adm/3/BGD_adm2.shp')
sf = os.path.dirname('/Users/ihasan/Downloads/Rivers of Bangladesh/natural/natural.shp')
Bangladesh = gpd.read_file(sdir1)
fr = gpd.read_file(sf)

Making the map of geospatial data

Now we will start to plot the map. As we already have loaded the necessary shapefiles. Now plotting is not a big deal. We will do some other stuff too. Such as how to plot a subset map, how to add colorbar, how to label shapefile. Ultimate product will look like this.

The final map showing all the geospatial data
The final map

Selecting partial data set

We will not plot all 64 districts. We will plot 7 specific districts and the natural forest in it. The district names are in a column of attribute table titled as ‘NAME_2’. The column name is saved in a variable named ‘cln’. The desired 7 districts names and the row number of that districts in the attribute table are saved in ‘std_a’ and ‘ids’ respectively.

cln = 'NAME_2'
std_a = ['Shatkhira', 'Khulna','Bagerhat', 'Jessore', 'Narail', 'Gopalgonj','Pirojpur']
ids = [np.where(Bangladesh[cln]==nm)[0] for nm in std_a]
ids = np.hstack(ids)

All seven districts of the main shapefiles are saved separately in variable name ‘studyar’. The plotting of that shapefile saved in a variable named ‘std_a1’. As it will be used later on.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)

The initial barebone map looks like this below.

Initial Map after visualising geospatial data
Initial Map

Plotting point shapefile

Now we will add the point shapefile that we have made earlier. The color of the point shapefile will vary according to the attribute table value of a column where a paremeter of each site was saved. The code adding point shapefile is given below.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)

After adding the point shapefile in the previous map the map will look like this below after mapping with python.

Map with geospatial data of polygon and point shapefile
Map with polygon and point shapefile

Resizing the mapping extant

After that we will add the forest shapefile in this region. The code and map is given below.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)
fr[fr['type']=='forest'].plot(ax=std_a1, color='g')
Map with forest
Map with forest

The forest shapefile has bigger extant than our study area. To resolve this we will limit x and y axis. The code and map will look like this below after limiting both axis.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)
fr[fr['type']=='forest'].plot(ax=std_a1, color='g')
plt.xlim(88.8, 90.4)
plt.ylim(21.5, 23.5)
Map with forest in it
The map with forest in it

Labelling the shapefile

Now its time for labelling the shapefiles. We will label the forest and the districts as well. For forest, one label is enough. But for districts, we need a for loop. We will also make white outline to all the label text.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)
fr[fr['type']=='forest'].plot(ax=std_a1, color='g')

tx1 = plt.text(x=89.4, y=22.25, s='Sundarban', fontsize=12)
tx1.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

for idx, row in studyar.iterrows():
    fnt = 10
    txt=plt.annotate(s=row[cln], xy=(studyar.geometry.centroid.x[idx],studyar.geometry.centroid.y[idx]),
                 horizontalalignment='center', fontsize=15, wrap=True, color='k')
    txt.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

plt.xlim(88.8, 90.4)
plt.ylim(21.5, 23.5)
Map with labels
Map with labels

Highlighting the area of interest

We will add a rectangle to specify our area of interest.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)
fr[fr['type']=='forest'].plot(ax=std_a1, color='g')

tx1 = plt.text(x=89.4, y=22.25, s='Sundarban', fontsize=12)
tx1.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

for idx, row in studyar.iterrows():
    fnt = 10
    txt=plt.annotate(s=row[cln], xy=(studyar.geometry.centroid.x[idx],studyar.geometry.centroid.y[idx]),
                 horizontalalignment='center', fontsize=15, wrap=True, color='k')
    txt.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

plt.plot([88.9, 90, 90, 88.9, 88.9],[21.6, 21.6, 23, 23,21.6], lw=4)
plt.xlim(88.8, 90.4)
plt.ylim(21.5, 23.5)
Map with rectangle
Map with rectangle

Adding colorbar

Now we will add the colour bar for the colour graduation of the point shapefile and a label to the colorbar.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)
fr[fr['type']=='forest'].plot(ax=std_a1, color='g')

tx1 = plt.text(x=89.4, y=22.25, s='Sundarban', fontsize=12)
tx1.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

for idx, row in studyar.iterrows():
    fnt = 10
    txt=plt.annotate(s=row[cln], xy=(studyar.geometry.centroid.x[idx],studyar.geometry.centroid.y[idx]),
                 horizontalalignment='center', fontsize=15, wrap=True, color='k')
    txt.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

plt.plot([88.9, 90, 90, 88.9, 88.9],[21.6, 21.6, 23, 23,21.6], lw=4)
plt.xlim(88.8, 90.4)
plt.ylim(21.5, 23.5)

sm = plt.cm.ScalarMappable(cmap='rainbow', norm=plt.Normalize(vmin=df['RE_all'].min(), vmax=50))
sm._A = []
cbar_axis = inset_axes(std_a1, width='10%', height='5%',loc='lower left',
                       bbox_to_anchor=[90.1, 22.2, 2, 0.8],bbox_transform=std_a1.transData)
cbar = plt.colorbar(sm, cax=cbar_axis, orientation='horizontal', pad=0.02)
std_a1.text(x=90.1, y=22.28, s='Recovery Efficiency \n of sites (%)', fontsize=9.8)
Map with colorbar and labels to it
Map with colorbar and labels to it

Subset map

Now we will add a subset map of whole Bangladesh to get an idea of our position of area of interest.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)
fr[fr['type']=='forest'].plot(ax=std_a1, color='g')

tx1 = plt.text(x=89.4, y=22.25, s='Sundarban', fontsize=12)
tx1.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

for idx, row in studyar.iterrows():
    fnt = 10
    txt=plt.annotate(s=row[cln], xy=(studyar.geometry.centroid.x[idx],studyar.geometry.centroid.y[idx]),
                 horizontalalignment='center', fontsize=15, wrap=True, color='k')
    txt.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

plt.plot([88.9, 90, 90, 88.9, 88.9],[21.6, 21.6, 23, 23,21.6], lw=4)
plt.xlim(88.8, 90.4)
plt.ylim(21.5, 23.5)


sm = plt.cm.ScalarMappable(cmap='rainbow', norm=plt.Normalize(vmin=df['RE_all'].min(), vmax=50))
sm._A = []
cbar_axis = inset_axes(std_a1, width='10%', height='5%',loc='lower left', bbox_to_anchor=[90.1, 22.2, 2, 0.8],bbox_transform=std_a1.transData)
cbar = plt.colorbar(sm, cax=cbar_axis, orientation='horizontal', pad=0.02)
std_a1.text(x=90.1, y=22.28, s='Recovery Efficiency \n of sites (%)', fontsize=9.8)

a=inset_axes(std_a1, width="30%", height='30%',loc='lower right')
Bangladesh.plot(figsize=(20,20), cmap='PuOr',ax=a)
a.text(x=88.6, y=22.25, s='Study Area', fontsize=9.8)
a.plot([88.9, 90, 90, 88.9, 88.9],[21.5, 21.5, 23, 23,21.5], lw=2)
Map with a subset map
Map with a subset map

Removing subset axis label and adding north arrow

For final touch, we will remove the axis tick label of the subset map, as this looks a bit mess and will add a north arrow to the map.

studyar = Bangladesh.iloc[ids,:]
std_a1=studyar.plot(figsize=(15,15),column=cln, cmap='Set3', edgecolor='k', lw=0.7)
csm = sites.plot(ax=std_a1, column='RE_all', cmap='rainbow', markersize=55)
fr[fr['type']=='forest'].plot(ax=std_a1, color='g')

tx1 = plt.text(x=89.4, y=22.25, s='Sundarban', fontsize=12)
tx1.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

for idx, row in studyar.iterrows():
    fnt = 10
    txt=plt.annotate(s=row[cln], xy=(studyar.geometry.centroid.x[idx],studyar.geometry.centroid.y[idx]),
                 horizontalalignment='center', fontsize=15, wrap=True, color='k')
    txt.set_path_effects([PathEffects.withStroke(linewidth=4, foreground='w')])

plt.plot([88.9, 90, 90, 88.9, 88.9],[21.6, 21.6, 23, 23,21.6], lw=4)

std_a1.text(x=90.214-0.025, y=22.55, s='N', fontsize=30)
std_a1.arrow(90.215, 22.36, 0, 0.18, length_includes_head=True,
          head_width=0.08, head_length=0.2, overhang=.1, facecolor='k')

plt.xlim(88.8, 90.4)
plt.ylim(21.5, 23.5)

sm = plt.cm.ScalarMappable(cmap='rainbow', norm=plt.Normalize(vmin=df['RE_all'].min(), vmax=50))
sm._A = []
cbar_axis = inset_axes(std_a1, width='10%', height='5%',loc='lower left',bbox_to_anchor=[90.1, 22.2, 2, 0.8], bbox_transform=std_a1.transData)
cbar = plt.colorbar(sm, cax=cbar_axis, orientation='horizontal', pad=0.02)
std_a1.text(x=90.1, y=22.28, s='Recovery Efficiency \n of sites (%)', fontsize=9.8)

a=inset_axes(std_a1, width="30%", height='30%',loc='lower right')
Bangladesh.plot(figsize=(20,20), cmap='PuOr',ax=a)
a.text(x=88.6, y=22.25, s='Study Area', fontsize=9.8)
a.plot([88.9, 90, 90, 88.9, 88.9],[21.5, 21.5, 23, 23,21.5], lw=2)
a.set_xticks([])
a.set_yticks([])
The final map
The final map

All the data and code regarding this post can be found in the GitHub repository link below.

https://github.com/hasanbdimran/Mapping-in-python

If some one need any help or tips or assistance regarding making map with python, you are always welcome to contact me through the comment section or via email of the contact page. If you want to use my expertise on making map or geoprocessing with python, you are most welcome to visit my fiverr gig about making map and geoprocessing.

https://www.fiverr.com/hasanbdimran/make-gis-map-using-python

The main objective of this blog or the fiverr gig is to explore my knowledge and also to learn from the expert visitors on this blog.

This Post Has 6 Comments

  1. Wasif

    Thanks for sharing. Keep it up and best wishes.

    1. admin

      Thanks for reading.

  2. Shishir

    Please Share the Upazilla Shape file

  3. Philipe Riskalla Leal

    Dear Mohammad Imran,

    I found your blog most helpful. I had never seen a more constructive example of the spatial plotting functions of python before.

    I was wondering, if perhaps you have a solution, or example, for inserting a scalebar to the plot. In most scientific institutions, the scalebar is mandatory, and it’s implementation is not trivial.

    Back in python V2.x, there was the option of the Basemap. Nevertheless, the migration to cartopy still lacks some functions previously available in Basemap.

    I thank you for your time, and I hope hearing from you soon.

    Sincerely,

    1. admin

      Thanks for your kind words. I have done some coding regarding on scale bar. But those are totally manual process. I will make those examples like general functions so that others can use and post on my blog soon.

Leave a Reply