[Python] Conversion from WGS84 to plane orthogonal coordinate system

Introduction

In surveying, the result is expressed not by latitude and longitude but by the coordinate values of a two-dimensional plane called the plane orthogonal coordinate system. It is a memo of how to perform the mutual conversion.

Currently, 19 Cartesian coordinate systems are used in Japan. ("Easy-to-understand plane orthogonal coordinate system" from Geospatial Information Authority of Japan)

Common (?) Misconception: Different from the coordinate system used in inertial navigation!

Coordinate system transformation sounds easy, but this is

--Convert the polar coordinates (latitude, longitude, ellipsoid height) of the ellipsoid of WGS84 to the coordinate values XYZ of the Cartesian coordinate system of the Earth center coordinate system (ECEF). --Rotate the values in the ECEF Cartesian coordinate system to correspond to the tangent plane in the corresponding ellipsoid.

It is "different" from the story. [Gauss-Krügel projection](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%E3%83%BB%E3%82%AF % E3% 83% AA% E3% 83% A5% E3% 83% BC% E3% 82% B2% E3% 83% AB% E5% 9B% B3% E6% B3% 95) However, I couldn't get there because the formula was complicated.

At first, I misunderstood it as a simple rotation calculation like the one above, and I was embarrassed to think back when I was talking to people around me saying, "It's such a simple calculation."

Support of Geographical Survey Institute

You can use the one provided by the Geospatial Information Authority of Japan for conversion. (As expected, a public institution!)

First, you can enter the coordinate values directly into the browser to convert.

https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xyf.html

It can also be used from the WEB API. Here are the features and descriptions provided. There is also a note that there is a limit of 10 times in 10 seconds from the same address so that access is not concentrated.

https://vldb.gsi.go.jp/sokuchi/surveycalc/main.html

This is an example of extracting by specifying the number of the plane orthogonal coordinate system.

use_gsi_webapi.py

import requests
from numpy import round
lat, lon = 35.71, 139.74

url = "http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xy.pl?"
params={'latitude':round(lat, 8), 'longitude': round(lon, 8), "refFrame":2, "zone": 9, 'outputType':'json'}

res =　requests.get(url, params=params)
if res.status_code == requests.codes.ok:
    print(res.json())

Here the following dictionary is returned:

{'OutputData': {
  'publicX': '-32170.0990',
  'publicY': '-8445.1361',
  'gridConv': '0.054477778',
  'scaleFactor': '0.99990088'}}

EPSG

The number of the plane orthogonal coordinate system in Japan is set by Japan, but it is set as EPSG internationally. Internationally, it is referred to by EPSG instead of the plane orthogonal coordinate system number. (For example, when using photogrammetry software metashape, etc.)

Although it is a source of information, I do not understand well the origin of the head family and the management organization. The following feels like authority.

-You can read the explanation at epsg.io. For example, WGS84. -Search for EPSG Geodetic Parameter Registry --The page spatial reference also has a list and search function.

Use pyproj

I tried to use the WEB API earlier, but this can also be done using the python library. I don't remember having a hard time installing it.

test_webapi.py

from pyproj import Transformer

lat, lon = 35.71, 139.74
wgs84_epsg, rect_epsg = 4326, 6677
tr = Transformer.from_proj(wgs84_epsg, rect_epsg)
x, y = tr.transform(lat, lon)
print(x, y)

When run, you get:

-32170.09900022997 -8445.136058616452

I don't know the significant figures, but they are in agreement with the results of the Geographical Survey Institute. (4 digits after the decimal point, 0.1 mm!)

Use of CRS

In the streets, the method of using CRS was introduced, but when I tried using it, ...

test_use_pyproj.py

import pyproj
wgs84=pyproj.CRS("EPSG:4326")
jgd2011_9 = pyproj.CRS("EPSG:6677")
lat, lon = 35.71, 139.74
print( pyproj.transform(wgs84, jgd2011_9, lat, lon) )

I get a warning when I run it, but I get the same result.

/usr/bin/ipython3:1: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
(-32170.09900022997, -8445.136058616452)

I think this is probably the old way. I have to read the tutorial carefully! https://pyproj4.github.io/pyproj/stable/examples.html

Summary / Miscellaneous feelings

--The conversion between WGS84 and the plane orthogonal coordinate system can be done with the API provided by the Geographical Survey Institute, but it seems that the same can be done with pyproj. ――Pyproj is multifunctional and has some interesting functions such as 4D tansform, so I would like to investigate if there is spare capacity.

Etc. (2020/08/08)