[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

mesh1.png

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

leaflet1.png

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)

leaflet2.png

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

static2.png

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

static3.png

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

static4.png 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))

static5.png 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

Visualize and understand Japan's regional mesh on a map
Visualize B League goals and goals on a heat map
Folium: Visualize data on a map with Python
Visualize grib2 on a map with python (matplotlib)
[Python] Visualize overseas Japanese soccer players on a map as of 2021.1.1
Visualize prefectures with many routes by prefecture on a Japanese map
A memo with Python2.7 and Python3 on CentOS
Map rent information on a map with python
Understand python lists, dictionaries, and so on.
[pyqtgraph] Understand SignalProxy and create a crosshair that follows the cursor on the graph
Create a web map using Python and GDAL
[Rails] How to display multiple markers on Google Map and display a balloon when clicked