[PYTHON] Visualize and understand Japan's regional mesh on a map
What is a regional mesh?
** Regional mesh ** is a mesh that has almost the same shape as the region based on latitude and longitude. One mesh is approximately square and is defined in various sizes. To summarize the things you use often,
-** Primary mesh : Approximately 80km on a side. - Secondary mesh : Approximately 10km on a side. The primary mesh is divided into 8 parts vertically and horizontally. - Tertiary mesh : Approximately 1km on a side. The secondary mesh is divided into 10 parts vertically and horizontally. - 4th mesh : Approximately 500m on a side. 3rd mesh divided into 2 parts vertically and horizontally. - 5th mesh **: Approximately 250m on a side. 4th mesh divided vertically and horizontally.
Reference: Statistics Bureau, Ministry of Internal Affairs and Communications | Regional Mesh Statistics
The code system is fixed for the mesh, and it is possible to calculate from latitude and longitude (Reference: Standard mesh system and code) .. However, the definition of each mesh code is also provided as data (e-Stat | Statistics on Map (Statistics GIS). -search? page = 1 & type = 2)). In this article, I would like to understand by using this to visualize the regional mesh.
Software library
Use Python (Anaconda / miniconda distribution). The required libraries can be installed from the following command.
conda install -c conda-forge jupyter geopandas descartes shapely \
matplotlib pip requests pillow chardet mplleaflet && \
pip install tilemapbase
The following libraries are characteristic of this purpose.
-geopandas: Used for processing geographic information data -mplleaflet: Plot on an interactive map -tilemapbase: Plot on a static map
Data acquisition
This time, we will use the data of the primary mesh 5339. The file can be downloaded directly from this link or Download the "M5339" shape file from the e-Stat page. I will. Since it is a Zip file, unzip it and if it has the following structure, it is ready.
└── QDDSWQ5339
├── MESH05339.dbf
├── MESH05339.prj
├── MESH05339.shp
└── MESH05339.shx
Data reading
Loading shapefiles is very easy with the geopandas library.
import geopandas as gpd
x = gpd.read_file("QDDSWQ5339/MESH05339.shp")
print(x.shape)
x.head()
(100800, 8)
KEY_CODE MESH1_ID MESH2_ID MESH3_ID MESH4_ID MESH5_ID OBJ_ID \
0 5339000011 5339 00 00 1 1 1
1 5339000012 5339 00 00 1 2 2
2 5339000013 5339 00 00 1 3 3
3 5339000014 5339 00 00 1 4 4
4 5339000021 5339 00 00 2 1 5
geometry
0 POLYGON ((139.00312 35.33333, 139.00000 35.333...
1 POLYGON ((139.00625 35.33333, 139.00312 35.333...
2 POLYGON ((139.00312 35.33542, 139.00000 35.335...
3 POLYGON ((139.00625 35.33542, 139.00312 35.335...
4 POLYGON ((139.00937 35.33333, 139.00625 35.333...
--The data contains all 5th order meshes starting with 5339.
--KEYCODE: The entire code of the mesh
--MESHX_ID: Code for the Xth mesh part. KEY_CODE is a combination of these from 1st to 5th.
--ʻOBJ_ID: The line number has no particular meaning (I think) --geometry`: Geographic information of the mesh. It is defined as a polygon.
Aggregation of polygons
Since the data this time is the data of the 5th order mesh, for example, the whole "5339" 1st order mesh itself is not defined anywhere. For coarse meshes, it can be obtained by aggregating 5th order polygons. The dissolve method of geopandas.GeoDataFrame is useful for this.
The following function aggregates the data in the 5th order mesh to the lower order data.
def aggregate_mesh(x, level):
tmp = x.copy()
code_len = [4, 6, 8, 9, 10][level-1]
tmp["key"] = tmp["KEY_CODE"].str.slice(0, code_len)
tmp = tmp[["key", "geometry"]]
tmp = tmp.dissolve(by="key")
return tmp
--code_len: Code length with specified particle size
--Specify the grain size code specified in tmp [" key "]
--Aggregate polygons with .dissolve
--. Dissolve takes a union for geographic information. I haven't used it this time, but if there are other columns, it will do calculations like groupby.aggregate at the same time.
# aggregate to mesh level 1
mesh1 = aggregate_mesh(x, 1)
mesh1
geometry
key
5339 POLYGON ((139.29062 35.33333, 139.28750 35.333...
This is the result of aggregating into the primary mesh. This time, only the primary mesh of "5339" is used, so it will be aggregated into one line.
Visualization of mesh geographic information (interactive map)
The geopandas.plotting module provides useful functions for visualizing geographic information. Especially for polygons, plot_polygon_collection is useful. Literally, it is a function that visualizes multiple polygons.
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
plot_polygon_collection(ax, mesh1.geometry, edgecolor="#121212", facecolor="#333333")
--Specify the drawing destination (ʻAxesofmatplotlib) with the first argument of plot_polygon_collection. So create ʻAxes from the pyplot.subplots function.
--The second argument is geopandas.GeoSeries
The result is a graph like the one above. I don't know where it is. Use the mplleaflet library to plot this on a map. This library provides the ability to draw matplotlib graphs on an interactive map.
import mplleaflet
mplleaflet.display(fig) #For jupyter
# mplleaflet.show(fig) #Other
When executed with Jupyter, the map will be displayed in the output as shown above. It can be enlarged and moved like a normal map app.
If you use .show, another page of the browser will open and a similar map will be displayed there.
You can see that the "5339" mesh is a mesh that includes a part of Tokyo to South Kanto. The reason why it is not exactly square is probably because the 5th order mesh of the sea part is not defined.
Similarly, let's display the secondary mesh.
mesh2 = aggregate_mesh(x, 2)
fig, ax = plt.subplots()
plot_polygon_collection(ax, mesh2.geometry, edgecolor="#121212", facecolor="#333333")
mplleaflet.display(fig)
Visualization of mesh geographic information (static map)
While mplleaflet is a very useful feature, it also has some weaknesses.
――The larger the amount to plot, the heavier it becomes. --At this time, it does not support drawing characters.
Regarding the first, for example, if you give a scatter plot of tens of thousands to 10,000, it will be difficult to draw.
As an alternative, try drawing on a static map. There seem to be some options, but I will use the TileMapBase library. This package takes an OpenStreetMap map as an image and draws it in the background. You can freely use the functions of matplotlib by overlaying the drawing on this.
The weakness is that the syntax is not as simple as mplleaflet, and the coordinate transformation must be done by the user.
It will be a little longer, but it is a function that draws a polygon on the map and draws a character string in the center.
import tilemapbase
tilemapbase.init(create=True)
from shapely import geometry
def plot_on_map_with_text(geoseries, texts):
# find the graph range
rect = geoseries.total_bounds
edgex = rect[2] - rect[0]
edgey = rect[3] - rect[1]
extent = tilemapbase.Extent.from_lonlat(
rect[0]-edgex*0.3, rect[2]+edgex*0.3,
rect[1]-edgey*0.3, rect[3]+edgey*0.3)
extent = extent.to_aspect(1.0)
fig, ax = plt.subplots(figsize=(8, 8), dpi=100)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
t = tilemapbase.tiles.build_OSM()
plotter = tilemapbase.Plotter(extent, t, width=600)
plotter.plot(ax, t)
polygons = []
centers = []
bounds = geoseries.bounds
for i in range(len(bounds)):
# convert to the plottable scale
minx, miny = tilemapbase.project(bounds['minx'][i], bounds['miny'][i])
maxx, maxy = tilemapbase.project(bounds['maxx'][i], bounds['maxy'][i])
polygons.append(
geometry.box(minx, miny, maxx, maxy))
centers.append([(minx + maxx)/2.0, (miny + maxy)/2.0])
polygons = gpd.GeoSeries(polygons)
plot_polygon_collection(ax, polygons,
edgecolor='#121212', facecolor='#000000', alpha=0.4)
for center, txt in zip(centers, texts):
ax.text(center[0], center[1], txt, fontdict={'color':'lightblue'})
return fig, ax
--tilemapbase.init needs to be executed once before use, and creates a cache of the acquired map image to avoid retrieving the same information repeatedly. It seems to be a consideration not to bother OpenStreetMap by repeating access.
--shapely is a library that provides the functions of basic elements of geographic information such as polygons, and geopandas also depends on it. Here, it is used to redefine the polygon whose coordinate system has been converted.
Plot the quadratic mesh with the last two digits of the code.
plot_on_map_with_text(mesh2.geometry, mesh2.index.str.slice(4, 6))
Looking at this, you can see the rules of mesh code. With the southwestern end as (0, 0), the first number increases when going north, and the second number increases when going east. Since the secondary mesh divides the primary mesh into eight, numbers from 0 to 7 appear.
You can draw in the same way after the 3rd order. However, the graph will be too detailed, so limit it to a specific area.
mesh3 = aggregate_mesh(x[x.MESH2_ID == "45"], 3)
plot_on_map_with_text(mesh3.geometry, mesh3.index.str.slice(6, 8))
The rules of the code haven't changed, and this time it's divided into 10 parts.
mesh4 = aggregate_mesh(x[(x.MESH2_ID == "45") & (x.MESH3_ID == "09")], 4)
plot_on_map_with_text(mesh4.geometry, mesh4.index.str.slice(8, 9))
The quaternary mesh divides the quaternary mesh into two parts vertically and horizontally, and assigns 1, 2, 3, and 4 codes to the southwest, southeast, northwest, and northeast, respectively.
mesh5 = aggregate_mesh(x[(x.MESH2_ID == "45") & (x.MESH3_ID == "09")], 5)
plot_on_map_with_text(mesh5.geometry, mesh5.index.str.slice(8, 10))
The four quaternary meshes are divided into five and displayed. The coding rules for the 5th mesh are the same as for the 4th mesh.
Recommended Posts