[Python] Conversion de WGS84 en système de coordonnées orthogonales plan

introduction

En topographie, le résultat n'est pas exprimé par la latitude et la longitude mais par les valeurs de coordonnées d'un plan bidimensionnel appelé système de coordonnées orthogonales du plan. Ceci est un mémo sur la façon d'effectuer la conversion mutuelle.

Actuellement, 19 systèmes de coordonnées orthogonales sont utilisés au Japon. ("Système de coordonnées orthogonales planes faciles à comprendre" de National Land Research Institute)

Malentendu courant (?): Différent du système de coordonnées utilisé dans la navigation inertielle!

La conversion du système de coordonnées semble facile, mais c'est

--Convertissez les coordonnées polaires (latitude, longitude, hauteur de l'ellipse) de l'ellipse de WGS84 en valeurs de coordonnées XYZ du système de coordonnées orthogonales du système de coordonnées du centre de la terre (ECEF).

• Faites pivoter les valeurs dans le système de coordonnées orthogonales ECEF pour qu'elles correspondent au plan tangent dans l'ellipse correspondante

C'est "différent" de l'histoire. [Projection Gauss-Kruger](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) Cependant, je n'ai pas pu y arriver car la formule était compliquée.

Au début, je l'ai mal compris comme un simple calcul de rotation comme celui ci-dessus, et j'ai été gêné de repenser quand je parlais aux gens autour de moi en disant: «C'est un calcul si simple».

Soutien du National Land Research Institute

Vous pouvez utiliser celui fourni par le National Land Institute pour la conversion. (Comme prévu, une institution publique!)

Tout d'abord, vous pouvez entrer les valeurs de coordonnées directement dans le navigateur pour la conversion.

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

Il peut également être utilisé à partir de l'API WEB. Voici les fonctionnalités et descriptions fournies. Il y a aussi une note qu'il y a une limite de 10 fois en 10 secondes à partir de la même adresse afin que l'accès ne soit pas concentré.

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

Voici un exemple d'extraction en spécifiant le numéro du système de coordonnées orthogonales plan.

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

Ici, le dictionnaire suivant est renvoyé:

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

EPSG

Les numéros de système de coordonnées orthogonales du plan du Japon sont définis par le Japon, mais au niveau international, ils sont définis comme EPSG. Au niveau international, il est désigné par EPSG au lieu du numéro du système de coordonnées rectangulaires plan. (Par exemple, lors de l'utilisation de metashape, qui est un logiciel d'arpentage photo.)

Bien que ce soit une source d'information, je ne comprends pas bien l'origine du chef de famille et de l'organisation de gestion. Ce qui suit ressemble à une autorité.

-Vous pouvez lire l'explication sur epsg.io. Par exemple, WGS84. -Vous pouvez rechercher par EPSG Geodetic Parameter Registry

• La page référence spatiale a également une fonction de liste et de recherche.

Utilisez pyproj

J'ai essayé d'utiliser l'API WEB plus tôt, mais cela peut également être fait en utilisant la bibliothèque python. Je ne me souviens pas avoir eu du mal à l'installer.

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)

Lors de l'exécution, vous obtenez:

-32170.09900022997 -8445.136058616452

Je ne connais pas les chiffres valides, mais ils sont en accord avec les résultats du Geographical Survey Institute. (4 chiffres après la virgule décimale, 0,1 mm!)

Utilisation de CRS

Dans les rues, comment utiliser CRS a été introduit, mais quand j'ai essayé de l'utiliser, ...

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

Je reçois un avertissement lorsque je l'exécute, mais j'obtiens le même résultat.

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

Je pense que c'est probablement l'ancienne méthode. Je dois lire attentivement le tutoriel! https://pyproj4.github.io/pyproj/stable/examples.html

Résumé / Sentiments divers

--La conversion entre WGS84 et le système de coordonnées orthogonales plan peut être effectuée avec l'API fournie par l'Institut géographique, mais il semble que la même chose puisse être faite avec pyproj. ――Pyproj est multifonctionnel et possède des fonctions intéressantes telles que le tansform 4D, je voudrais donc vérifier s'il y a une capacité de réserve.

Etc. (08/08/2020)