Square mesh coordinates of Japan
Introduction
We may want to know the coordinates at a fixed distance over the entire land of Japan because we want to find the land coverage of wireless communication systems, or to investigate the radiowave propagation characteristics that depend on the region.
As for the latitude, which is the north-south angle from the center of the earth to the equator, one circumference of the earth is about 40,000 kilometers, so the distance per degree is about 111 kilometers. On the other hand, regarding the longitude, which is the angle in the east-west direction from the Greenwich Observatory in the United Kingdom, the distance per degree depends on the latitude, and it is the value obtained by multiplying the previous 111 kilometers by the cosine of the latitude.
So I tried to draw a line for each latitude and longitude on a large and easy-to-coordinate map such as the segmented air route map, and to read the entire coordinate range of the land. However, I soon realized that this needs a lot of work.
As I was thinking about alternatives, I realized that I could use the regional mesh code of the Statistics Bureau of MIC (the Ministry of Internal Affairs and Communications).
Regional mesh code
A region mesh is a mesh of regions of approximately the same size, based on latitude and longitude. The mesh code published by the Statistics Bureau is in CSV (comma separated values) format for each region of 47 prefectures (however, Hokkaido is so large and is divided into three).
The mesh code is quantifies by the latitude and longitude and hierarchized. The order of the mesh code is also defined according to the size of one side. In this mesh code, the primary mesh is shown in 4 digits with a side of 80 kilometers, the secondary mesh is shown in 6 digits with a side of 10 kilometers, and the tertiary mesh is shown in 8 digits with a side of 1 kilometer. The tertiary mesh code excluding the last two digits is equivalent to the secondary mesh code. The area mesh code includes the “prefecture, city, town, and village code of Japan” (5 digits), as well as the tertiary mesh code. This mesh code includes areas that are extremely difficult to visit even for Japanese people.
First, I download all of this regional mesh CSV file. The formula for converting the regional mesh code to coordinates is written on page 13 of the characteristics and history of regional mesh statistics of Japan (PDF: 1,249KB). I used the Python code shown in the article and I made it as a function.
mesh2latlon.py
#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# https://www.gis-py.com/entry/py-mesh
import sys
def mesh2latlon (meshCode):
code_first_two = meshCode[0:2] # for the primary mesh code
code_last_two = meshCode[2:4]
code_first_two = int(code_first_two)
code_last_two = int(code_last_two)
lat = code_first_two * 2 / 3
lon = code_last_two + 100
if len(meshCode) > 4: # for the secondary mesh code
if len(meshCode) >= 6:
code_fifth = meshCode[4:5]
code_sixth = meshCode[5:6]
code_fifth = int(code_fifth)
code_sixth = int(code_sixth)
lat += code_fifth * 2 / 3 / 8
lon += code_sixth / 8
if len(meshCode) == 8: # for the tertiary mesh code
code_seventh = meshCode[6:7]
code_eighth = meshCode[7:8]
code_seventh = int(code_seventh)
code_eighth = int(code_eighth)
lat += code_seventh * 2 / 3 / 8 / 10
lon += code_eighth / 8 / 10
return lat, lon
if __name__ == '__main__':
if len (sys.argv) < 2:
print ("Usage {} meshCode".format(sys.argv[0]))
meshCode = '52396594'
lat, lon = mesh2latlon (meshCode)
print ("meshCode {} is converted to {:.5f} {:.5f}."
.format(meshCode, lat, lon))
sys.exit (1)
lat, lon = mesh2latlon(sys.argv[1])
print (lat, lon)
# EOF
Since the mesh required this time is 10 kilometers on a side, I create the Python code japanloc.py that outputs the coordinates while reading all the CSV files and that imports the code mesh2latlon/py.
japanloc.py
#! /usr/bin/env python3
# -*- coding: utf-8 -*-
import csv
import glob
import mesh2latlon
import os
import sys
import time
import unicodedata
fnames = glob.glob ("csv/*.csv")
fnames.sort()
for fname in fnames:
with open (fname, 'r', encoding='CP932') as fp:
line = csv.reader (fp)
header = next (line)
mcode_prev = ''
for row in line:
mcode = row[2]
if mcode[0:6] == mcode_prev[0:6]: continue # for use with the secondary mesh code
mcode_prev = mcode
loc = row[1]
try:
if row[3]: loc += row[3]
except IndexError:
pass
lat, lon = mesh2latlon.mesh2latlon (mcode)
# Japanese character width https://mulberrytassel.com/python-practice-zenhan/
# cnt = 0
# for c in loc:
# if unicodedata.east_asian_width (c) in "FWA": cnt += 1
# col = 28 - cnt
# print ("{:{wd}s} {} {:.2f} {:.2f}"
# .format (loc, mcode, lat, lon, wd = col))
print ("{} {:.6f} {:.6f}".format (mcode, lat, lon))
# EOF
The code commented is for displaying the city name and was used for debugging. This code outputs the secondary mesh code, latitude, and longitude line by line. The same mesh code is given to the boundaries of cities, towns and villages. To remove this duplicate mesh code, hide the city names, sort them in ascending order of mesh code, and use the uniq command to remove the duplicates.
./japanloc.py | sort -k 1n | uniq > japanloc.txt
These coordinates sorted by mesh code point to places that are far apart for each row, so I sort the column of the latitude as the first priority and the column of the longitude as the second priority. To do this sorting, I load it into vim to interactively check the result with a vim command :%! sort -k 2n -k 3n and save the result with the ZZ command.
The number of lines in the 10-kilometer square coordinate file japanloc.txt obtained in this way is 11,986 lines. The total number of lines in the tertiary mesh excluding the header in the first line of the Statistics Bureau CSV file was 461,373 lines, and the number of lines after extracting the secondary mesh was 12,879 lines.
Application of the mesh code
In fact, the reason for thinking about this is the coordinate map of the reference station of the ALES centimeter-level positioning service analyzed by Professor Tomoji Takasu of Tokyo University of Maritime Science and Technology. Then, I thought it might be a great tool of analyzing the RTK (reltime kinematic) TTFF (time to first fix) dependence of not only for a distance but also for an azimuth. I aware the service is useful not only for positioning but also for research.
Therefore, in order to know the coordinates of the reference station, I created a code to display RTCM (Radio Technical Commission for Maritime Services) message types 1005 (coordinates) and 1033 (station number, antenna name, antenna setting) the number, receiver name and serial number) on one line, while sending the above coordinate sequence as NMEA GGA to this NTRIP Caster. This result is recorded only when there is a change in the RTCM message. This code was created by importing showrtcm.py of the QZS L6 Tools and was in private.
Fortunately, the service transmits not only the observation data per second, but also transmits the coordinate values per second. After sending one our coordinate, I waited for 1.5 seconds and repeated the procedure to send the next coordinate, therefore it took approximately 5 hours. This result listed duplicate reference stations. As in the previous step, I removed the duplicate reference station entries with the uniq command after sorting in ascending order of station numbers.
By this method, I found 3,086 reference stations. Professor Takasu has discovered 3,194 reference stations. Therefore, my results seem to have missed more than 100 stations. For example, station number 172 was duplicated in two separate reference stations, and some reference stations were given the receiver version number 4.4.0.
Conclusion
The 10-kilometer square coordinates of Japan were obtained using the CSV file of the Statistics Bureau of MIC, Japan. Densely installed RTK reference stations can be a useful tool for positioning and navigation research. Since Japan has the quasi-zenith satellites and excellent commercial services, it is highly valuable to conduct positioning research and learn positioning technology in Japan.