[PYTHON] Modèle d'élévation numérique utilisant les informations cartographiques de base

introduction

Décrit comment utiliser le modèle numérique d'élévation (DEM) publié par le National Land Research Institute. Décrit comment convertir des données DEM en données d'élévation 2D.

nom de fichier

Données cartographiques de base (modèle d'élévation numérique) Données de maillage de 5 m (altitude) Le nom de fichier de est au format suivant.

FG-GML-pppp-qq-rr-DEM5X-yyyymmdd.xml

pppp est le numéro de maillage principal qq est le numéro de maillage secondaire rr est le numéro de maillage tertiaire X est une méthode d'enquête, A: levé laser, B: levé photo. aaaammjj est la date.

--Exemple

FG-GML-5338-06-89-DEM5A-20130702.xml

Numéro de maillage

Le maillage principal est un carré d'une largeur de 1 degré de longitude et d'une hauteur de 40 minutes de latitude. Le numéro de maillage principal est représenté par la latitude et la longitude de l'extrémité sud-ouest de ce carré. Le numéro de maillage principal est représenté par 4 chiffres. Les deux premiers chiffres correspondent à 1,5 fois la latitude. Les deux derniers chiffres sont les deux derniers chiffres de la longitude.

Le maillage secondaire est le maillage principal divisé verticalement et horizontalement en huit parties égales. Numéro (nord-sud, est-ouest) avec l'extrémité sud-ouest comme (0,0). La plage de numéros est comprise entre 0 et 7. Le numéro de maillage secondaire est un nombre à deux chiffres disposé dans l'ordre nord-sud et est-ouest.

Le maillage tertiaire est le maillage secondaire divisé verticalement et horizontalement en 10 parties égales. La numérotation est la même idée que le maillage secondaire.

image

Contenu du fichier

Seules les informations principales sont extraites

Il s'agit d'un tableau de numéros de maillage primaire, secondaire et tertiaire.

  <mesh>53380689</mesh>

Décrit le nombre de cellules de grille qui composent le maillage. L'attribut bas est un nombre indiquant l'extrémité nord-ouest. C'est toujours (0,0). L'attribut haut est un nombre indiquant l'extrémité sud-est. La maille de 5 m est fixée à (224,149) et la maille de 10 m est fixée à (1124,749).

  <gml:limits>
  	<gml:GridEnvelope>
  		<gml:low>0 0</gml:low>
  		<gml:high>224 149</gml:high>
  	</gml:GridEnvelope>
  </gml:limits>

Les données d'élévation de chaque cellule de la grille sont répertoriées. Les données incluses dans la balise tupleList sont les données d'élévation. S'il n'y a pas de données sur la surface de l'eau ou l'altitude, la valeur d'altitude sera "-9999.".

  <gml:DataBlock>
  	<gml:rangeParameters>

<gml: QuantityList uom = "DEM constituants points" /> </gml:rangeParameters> gml:tupleList (Données d'altitude) Exemple) Surface du sol, 1055,33 </gml:tupleList> </gml:DataBlock>

  			<gml:sequenceRule  order="+x-y">Linear</gml:sequenceRule>

Spécifiez la séquence de données avec sequenceRule. «+ X» dans l'ordre indique que x est dans la direction positive (ouest → est), et «-y» indique que y est dans la direction négative (nord → sud). De l'extrémité nord-ouest à l'extrémité sud-est, les données disposées en ligne comme indiqué ci-dessous sont placées dans la balise tupleList. 1.JPG

  <gml:startPoint>173 0</gml:startPoint>

Spécifiez le point de départ de la cellule de la grille avec startPoint. ** Les données placées dans le DataBlock seront les numéros de cellule de la grille (173,0) à (224,149). ** ** ** Veuillez noter que le nombre réel de données n'est pas nécessairement (gml: low + 1) * (gml: high + 1). ** **

Récupération de données

Maintenant que nous connaissons le contenu du XML, nous allons extraire les données DEM sous forme de données d'élévation bidimensionnelles. J'utiliserai Python pour le récupérer. Il s'agit d'un type ndarray 2D de Numpy. La taille du tableau est fixée au nombre de colonnes (gml: low + 1) et au nombre de lignes (gml: high + 1). Toutes les données avant startPoint sont définies sur «-9999».

hyoukou.py


import sys
import re
import numpy as np
import matplotlib.pyplot as plt

with open(sys.argv[1], "r") as f:
    # mesh
    r = re.compile("<mesh>(.+)</mesh>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            mesh = m.group(1)
            break
    # grid len
    r = re.compile("<gml:high>(.+) (.+)</gml:high>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            xlen = int(m.group(1)) + 1
            ylen = int(m.group(2)) + 1
            break
    # start
    r = re.compile("<gml:tupleList>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            break
    # data
    data = []
    r  = re.compile("</gml:tupleList>")
    r2 = re.compile(",(.+)")
    for ln in f:
        m = r.search(ln)
        if m != None:
            break
        else:
            m = r2.search(ln)
            if m != None:
                data.append(float(m.group(1)))
    # start point
    startx = starty = 0
    r = re.compile("<gml:startPoint>(.+) (.+)</gml:startPoint>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            startx = int(m.group(1))
            starty = int(m.group(2))
            break

data2 = np.empty(xlen*ylen)
start_pos = starty*xlen + startx

for i in range(xlen*ylen):
    if i < start_pos:
        data2[i] = -9999.
    else:
        data2[i] = data[i-start_pos]

data = data2.reshape(ylen, xlen)

#print(data)
plt.imshow(data)
plt.colorbar()
plt.show()

Comment appeler hyoukou.py Appelez avec le nom de fichier DEM comme argument. Affichez les données d'élévation 2D sous forme d'image à l'aide de matplotlib. S'il n'y a pas de perte de données (-9999.), Il peut être affiché normalement.

>python hyoukou.py FG-GML-5338-06-89-DEM5B-20130702.xml

Recommended Posts

Modèle d'élévation numérique utilisant les informations cartographiques de base
Informations sur la carte de base à l'aide de la conversion Geotiff Python des données numériques d'élévation