[PYTHON] Visualize location information with Basemap

There are many ways to visualize geospatial information using Python, but this time I will use ** Basemap Matplotlib Toolkit ** to visualize location information.

Click here for articles introduced in the past regarding geospatial information visualization, for example.

1. About Basemap

Simply put, Basemap allows you to add matplotlib plotting capabilities while drawing various map projections, map tiles, coastlines, rivers, and more. It is mainly used by geoscientists.

The matplotlib basemap toolkit is a library for plotting 2D data on maps in Python. Basemap does not do any plotting on it’s own, but provides the facilities to transform coordinates to one of 25 different map projections (using the PROJ.4 C library). Matplotlib is then used to plot contours, images, vectors, lines or points in the transformed coordinates. Shoreline, river and political boundary datasets (from Generic Mapping Tools) are provided, along with methods for plotting them. The GEOS library is used internally to clip the coastline and polticial boundary features to the desired map projection region. Basemap is geared toward the needs of earth scientists, particularly oceanographers and meteorologists. (Quoted from https://matplotlib.org/basemap/users/intro.html)

1-1. Basics of Basemap

First, use conda to install baseline.

$ conda install -c anaconda basemap
$ conda install -c conda-forge basemap-data-hires
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

fig = plt.figure()
#Specify the range to draw the map.
m = Basemap(llcrnrlat=30, urcrnrlat=50, llcrnrlon=125, urcrnrlon=150) #Create Basemap instance
#Subtract latitude / longitude every 10 degrees. The second option sets where to put the label up, down, left and right.
m.drawparallels(np.arange(-90, 90, 10), labels=[ True,False, True, False]) #Latitude line
m.drawmeridians(np.arange(0, 360, 10), labels=[True, False, False, True]) #Longitude line
m.drawcoastlines() #Border line
plt.show()

image.png

1-2. Detailed map settings

You can flexibly change the map resolution and projection method. (There are many detailed settings ...)

class mpl_toolkits.basemap.Basemap(llcrnrlon=None, llcrnrlat=None, urcrnrlon=None, urcrnrlat=None, llcrnrx=None, llcrnry=None, urcrnrx=None, urcrnry=None, width=None, height=None, projection='cyl', resolution='c', area_thresh=None, rsphere=6370997.0, ellps=None, lat_ts=None, lat_1=None, lat_2=None, lat_0=None, lon_0=None, lon_1=None, lon_2=None, o_lon_p=None, o_lat_p=None, k_0=None, no_rot=False, suppress_ticks=True, satellite_height=35786000, boundinglat=None, fix_aspect=True, anchor='C', celestial=False, round=False, epsg=None, ax=None)

projection\resolution Low resolution Medium resolution High resolution
Equirectangular projection image.png image.png image.png
Mercator projection image.png image.png image.png
Lambert projection image.png image.png image.png

By the way, the default value when nothing is set is projection ='cyl' (equirectangular projection) resolution ='c' (coarse resolution)

2. Location information visualization method

2-1. Data source

"Pseudo-human flow data" based on SNS analysis data by the Center for Spatial Information Science, University of Tokyo

This time, I would like to use the pseudo-human flow data that is released for free as a sample.

2-2. Data structure

Read the Kansai area pseudo-human flow data on September 16, 2013 and check the contents.

import pandas as pd
df = pd.read_csv('./Kansai/2013_09_16.csv')
df

'''
First column: User ID
Second column: Gender estimate (1):Male, 2:Female, 0.3: Unknown, NA: Unestimated)
3rd column: Date / time (24 hours every 5 minutes)
Fourth column: Latitude
5th column: Longitude
6th column: Resident category(Major classification)* Character string type
7th column: Resident category(Subcategory)* Character string type
8th column: State(Stay or move)* Character string type
9th column: Resident category ID(Corresponds to the 6th and 7th lines)
'''

This time, we will target the moving user (STAY_MOVE =='MOVE'). Also, while dividing the timestamp into hours and minutes, give a rank for each user and each hour and format it into a form that is easy to handle.

from dfply import *

df = df >> filter_by(X.STAY_MOVE=='MOVE') >> select(columns_to(X.lon, inclusive=True))
df = df >> separate(X.timestamp, ['col1','hour','col2','minute','col3'], sep=[10,13,14,16],convert=True) >> select(~X.col1, ~X.col2, ~X.col3)
df = df >> group_by(X.uid,X.hour) >> mutate(rk=row_number(X.minute))
df

image.png

2-3. Basemap initial value setting and location information visualization

Now, I'm going to visualize the location information, but first I will set the map that I will consider as the base this time. There are also methods that can draw the background map of ArcGIS and the borders of prefectures and cities, towns and villages, so I will apply them as well.

For prefectural borders and municipal borders, please download the Shapefile from __ here __. The file structure is as follows.

-gadm
 -gadm36_JPN_1.cpg 
 -gadm36_JPN_1.shp
 -gadm36_JPN_2.dbf
 -gadm36_JPN_2.shx
 -gadm36_JPN_1.dbf
 -gadm36_JPN_1.shx
 -gadm36_JPN_2.prj
 -gadm36_JPN_1.prj
 -gadm36_JPN_2.cpg
 -gadm36_JPN_2.shp


** Initial map setting **

def basemap(): 
    fig = plt.figure(dpi=300)
    m = Basemap(projection="cyl", resolution="i", llcrnrlat=33.5,urcrnrlat=36, llcrnrlon=134, urcrnrlon=137)
    m.drawparallels(np.arange(-90, 90, 0.5), labels=[True, False, True, False],linewidth=0.0, fontsize=8)
    m.drawmeridians(np.arange(0, 360, 0.5), labels=[True, False, False, True],linewidth=0.0, rotation=45, fontsize=8)
    m.drawcoastlines(color='k')
    m.readshapefile('./gadm/gadm36_JPN_1', 'prefectural_bound1', color='k', linewidth=.8) #Prefectural border
    m.readshapefile('./gadm/gadm36_JPN_2', 'prefectural_bound2', color='lightgray', linewidth=.5) #Municipal border
    m.arcgisimage(service='World_Street_Map', verbose=True, xpixels=1000, dpi=300)

image.png

** If you do so far, you can plot like matplotlib. ** **

https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#scatter

tmp1=df[(df['gender']==1) & (df['hour']==9) & (df['rk']==1)]
tmp2=df[(df['gender']==2) & (df['hour']==9) & (df['rk']==1)]

basemap()
plt.scatter(tmp1['lon'],tmp1['lat'],color='b',s=0.5) #Men in blue
plt.scatter(tmp2['lon'],tmp2['lat'],color='r',s=0.5) #Woman in red

Here, only the first log at 9 o'clock is plotted for each user. image.png

https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#hexbin

tmp=df[(df['hour']==9) & (df['rk']==1)]

basemap()
plt.hexbin(tmp['lon'],tmp['lat'], gridsize=50, cmap='rainbow', mincnt=1, bins='log',linewidths=0.3, edgecolors='k')

You can also change the size of the hexagrid with gridsize. image.png

3. Bonus

For example, it is possible to multiply weather data that is likely to affect the movement of people. By the way, September 16, 2013 was the day when the typhoon passed.

image.png out.gif
From the Japan Meteorological Agency 1km mesh radar echo

Vector display is also possible by calculating the speed and azimuth from the position information of the start point and end point. (Aside from whether there is a use) out.jpeg

That's it! Until the end Thank you for reading!

Recommended Posts

Visualize location information with Basemap
Visualize latitude / longitude coordinate information with kepler.gl
Quickly visualize with Pandas
Visualize data with Streamlit
Visualize claims with AI
termux × AWS Send smartphone location information to AWS with IoT
Visualize 2ch threads with WordCloud-Scraping-
Get information with zabbix api
Visualize Wikidata knowledge with Neo4j
Get Alembic information with Python
Get video file information with ffmpeg-python
Visualize decision trees with jupyter notebook
Visualize python package dependencies with graphviz
[python] Read information with Redmine API
Get weather information with Python & scraping
Notify LINE of location information (Google Maps) with GPS multi-unit SORACOM Edition