Transformez plusieurs données numériques d'élévation en une seule image avec Python
Suite de l'article que j'ai écrit plus tôt? Quelque chose comme
"Essayez d'imaginer les données d'altitude du National Land Research Institute avec Python"
Dans cet article, j'aimerais lancer plusieurs fichiers xml et les organiser correctement pour créer une seule image.
Préparation
Cette fois, les données de maillage de n'importe quelle préfecture seront entièrement supprimées.
Sélectionnez et achetez n'importe quel maillage à partir de la page de téléchargement du National Land Research Institute.
https://fgd.gsi.go.jp/download/menu.php
Cliquez sur le bouton "Aller à la sélection de fichier" dans "Informations cartographiques de base Modèle d'élévation numérique" pour passer à l'écran de téléchargement.
Cette fois, j'utiliserai les données de la préfecture d'Osaka.
Veuillez cocher 10 (ligne de contour de la carte topographique) de mailles de 10 m sur le côté gauche.
Dans le code écrit dans cet article, les données fournies par le National Land Research Institute sont en échelle de gris telles quelles
En d'autres termes, dans le cas de données de maillage de 10 m, une image en niveaux de gris de 1125 x 750 est générée pour chaque fichier xml, et ils sont disposés un par un pour générer une grande image.
Dans le cas de la préfecture d'Osaka, une grande image stupide de 6x1125px en largeur et 10x750px en hauteur sera générée, il est donc préférable de sélectionner une petite préfecture.
Si vous souhaitez réduire le nombre de pixels dans une seule image, veuillez jouer avec les parties décrites ci-dessous.
En passant, une fois que vous avez répondu au fichier ZIP déposé, je pense qu'il a la structure de données suivante.
Répertoire décompressé |-- FG-GML-0000-00-DEM10B.zip |-- FG-GML-0000-00-DEM10B.zip |-- etc...
Vous pouvez décompresser chaque fichier ZIP davantage, mais cette fois, il sera traité comme un fichier ZIP car c'est un peu gênant.
Le code entier ressemble à ci-dessous
import os
import re
import sys
import glob
import zipfile
import numpy as np
from PIL import Image
import xml.etree.ElementTree as ET
WIDTH = 1125
HEIGHT = 750
#Préparation
def set_value(dir, rate):
global WIDTH
global HEIGHT
#Répertorier et trier les maillages par nom de fichier
meshList = []
for f in glob.glob(os.path.join(dir, "*.zip")):
base = os.path.basename(f)
meshList.append(base[7:11] + base[12:14])
meshList.sort()
#Nombre de données à traiter
denominator = len(meshList)
#Nombre de données traitées
numerator = 0
#Examinez le maillage au nord, au sud, à l'est et à l'ouest dans la liste des mailles
top = right = meshList[-1]
bottom = left = meshList[0]
for meshNumber in meshList:
str(meshNumber)
if top[0:2] < meshNumber[0:2] or (top[0:2] <= meshNumber[0:2] and top[4] <= meshNumber[4] and top[4] <= meshNumber[4]):
top = meshNumber
if bottom[0:2] > meshNumber[0:2] or (bottom[0:2] >= meshNumber[0:2] and bottom[4] >= meshNumber[4]):
bottom = meshNumber
if right[2:4] < meshNumber[2:4] or (right[2:4] <= meshNumber[2:4] and right[5] <= meshNumber[5]):
right = meshNumber
if left[2:4] > meshNumber[2:4] or (left[2:4] >= meshNumber[2:4] and left[5] >= meshNumber[5]):
left = meshNumber
#Obtenez la taille de la toile à partir du maillage final
vartical = (int(top[0:2])-int(bottom[0:2])) * 8 + int(top[4]) - int(bottom[4]) + 1
horizon = (int(right[2:4])-int(left[2:4])) * 8 + int(right[5]) - int(left[5]) + 1
maxArea = vartical * horizon
point = top[0:2] + left[2:4] + top[4] + left[5]
#génération de canevas
canvas = Image.new('RGB', (int(WIDTH / rate) * horizon, int(HEIGHT / rate) * vartical), (0, 0, 0))
#Lecture de fichiers et imagerie
for zipfile in glob.glob(os.path.join(dir, "*.zip")):
paste_image(zipfile, point, canvas, rate)
numerator += 1
print('%d / %d' % (numerator, denominator))
canvas.save('dem.png', 'PNG', quality=100)
#Méthode d'imagerie
def paste_image(z, p, c, r):
global WIDTH
global HEIGHT
point = p
canvas = c
rate = r
with zipfile.ZipFile(z, "r") as zf:
for info in zf.infolist():
inner = info
with zf.open(inner) as zfile:
root = ET.fromstring(zfile.read())
namespace = {
'xml': 'http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema',
'gml': 'http://www.opengis.net/gml/3.2'
}
dem = root.find('xml:DEM', namespace)
#Numéro de maillage
mesh = dem.find('xml:mesh', namespace).text
#Nombre de cellules disposées(La valeur réelle est ajoutée par 1)
high = dem.find('./xml:coverage/gml:gridDomain/gml:Grid/gml:limits/gml:GridEnvelope/gml:high', namespace).text.split(' ')
highX = int(high[0]) + 1
highY = int(high[1]) + 1
#Réglage de la taille de l'image(Le nombre de données==Nombre de pixels)
dataSize = highX * highY
#Tableau de données d'altitude
dem_text = dem.find('./xml:coverage/gml:rangeSet/gml:DataBlock/gml:tupleList', namespace).text
data = re.findall(',(.*)\n', dem_text)
dataNp = np.empty(dataSize)
for i in range(len(data)):
if(data[i] == "-9999.00"):
dataNp[i] = 0
else:
dataNp[i] = float(data[i])
#Coordonnées de début des données
startPoint = dem.find('./xml:coverage/gml:coverageFunction/gml:GridFunction/gml:startPoint', namespace).text.split(' ')
startPointX = int(startPoint[0])
startPointY = int(startPoint[1])
startPosition = startPointY * highX + startPointX
#Lorsque le nombre de données est insuffisant(S'il y a des espaces au-dessus et en dessous de l'image)
if(len(dataNp) < dataSize):
add = []
#Lorsque les données en bas de l'image sont insuffisantes
if(startPosition == 0):
for i in range(dataSize - len(dataNp)):
add.append(0)
dataNp.extend(add)
#Lorsque les données en haut et en bas de l'image sont insuffisantes
else:
for i in range(startPosition):
add.append(0)
dataNp[0:0] = add
add = []
for i in range(dataSize - len(dataNp) - len(add)):
add.append(0)
dataNp.extend(add)
#Conversion d'entiers 8 bits de données
dataNp = (dataNp / 15).astype(np.uint8) #Divisez dataNp par 15 pour correspondre au point le plus élevé du mont Fuji à 255
data = dataNp.reshape(highY, highX)
#Imagerie des données numériques d'élévation
pilImg = Image.fromarray(np.uint8(data))
pilImg = pilImg.resize((int(highX / rate), int(highY / rate)), Image.LANCZOS) # NEAREST
#Calculer les coordonnées du maillage à coller
targetX = (int(point[0:2]) - int(mesh[0:2])) * 8 + int(point[4]) - int(mesh[4])
targetY = (int(mesh[2:4]) - int(point[2:4])) * 8 + int(mesh[5]) - int(point[5])
canvas.paste(pilImg, (targetY * int(WIDTH / rate), targetX * int(HEIGHT / rate)))
#Entrez le répertoire contenant le fichier ZIP
DIR = sys.argv[1]
RATE = int(sys.argv[2])
set_value(DIR, RATE)
Lors de l'exécution de python, cela devrait fonctionner si vous spécifiez un répertoire plein de fichiers ZIP
Supplément
Où générer un canevas avec la méthode set_value
canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))
Les 1125 et 750 de Coco écrivent chacun simplement manuellement le nombre de cellules dans les données numériques d'élévation.
Étant donné que les valeurs elles-mêmes sont les mêmes que les highX et highY qui apparaissent dans le code, c'est un tracas deux fois, et les highX et Y sont déjà divers codes merdiques tels que la recherche de la valeur chaque fois qu'une image est générée.
(Je vais le réécrire peut-être)
Je suis content que quelqu'un le refactore ...
Et
dataNp = (dataNp / 15).astype(np.uint8)
Qu'est-ce que 15? Veuillez vous référer à l'article précédent pour la réponse.
https://qiita.com/artistan/items/99266407702d4152e9d5
Redimensionnons la grande image
Comme je l'ai écrit au début, il n'est pas pratique de générer une grande image avec le code tel qu'il est.
Donc, la partie suivante de la méthode set_value
canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))
De
canvas = Image.new('RGB', (1125/10 * horizon, 750/10 * vartical), (0, 0, 0))
Et redimensionner le pilImage commenté
pilImg = pilImg.resize((int(highX)/ 10, int(highY) / 10), Image.LANCZOS) # NEAREST
Si vous le faites, il sera pratique de générer une image réduite à 1/10.
Enfin
Voici la préfecture d'Osaka, qui a été générée avec le code ci-dessus et réduite de manière appropriée.
Osaka n'est pas intéressante car elle n'a pas beaucoup de hauteur, mais je pense que cela vaut la peine d'être vu dans un endroit plus accidenté.
Sur une échelle de gris, la forme de la montagne ressemble à une nervure de feuille, qui est plutôt fraîche.
Je ne pouvais que l'expliquer incompréhensible, mais merci d'avoir lu.
Recommended Posts