Machen Sie mit Python mehrere numerische Höhendaten zu einem einzigen Bild
Fortsetzung des Artikels, den ich früher geschrieben habe? So etwas wie
"Versuchen Sie, die Höhendaten des National Land Research Institute mit Python abzubilden"
In diesem Artikel möchte ich mehrere XML-Dateien einfügen und sie richtig anordnen, um ein einzelnes Bild zu erstellen.
Vorbereitung
Dieses Mal werden die Maschendaten einer Präfektur vollständig gelöscht.
Wählen und erwerben Sie ein beliebiges Netz auf der Download-Seite des National Land Research Institute.
https://fgd.gsi.go.jp/download/menu.php
Klicken Sie unter "Grundlegende Karteninformationen Numerisches Höhenmodell" auf die Schaltfläche "Zur Dateiauswahl wechseln", um zum Download-Bildschirm zu springen.
Dieses Mal werde ich die Daten aus der Präfektur Osaka verwenden.
Bitte überprüfen Sie 10 (Konturlinie der topografischen Karte) von 10 m Maschenweite auf der linken Seite.
In dem in diesem Artikel geschriebenen Code sind die vom National Land Research Institute bereitgestellten Daten so wie sie sind grau skaliert
Mit anderen Worten, im Fall von 10-m-Maschendaten wird für jede XML-Datei ein Graustufenbild von 1125 x 750 erzeugt, und sie werden einzeln angeordnet, um ein großes Bild zu erzeugen.
Im Fall der Präfektur Osaka wird ein dummes großes Bild mit einer Breite von 6 x 1125 Pixel und einer Höhe von 10 x 750 Pixel erzeugt. Wählen Sie daher am besten eine kleine Präfektur aus.
Wenn Sie die Anzahl der Pixel in einem einzelnen Bild reduzieren möchten, spielen Sie bitte mit den unten beschriebenen Teilen.
Übrigens, sobald Sie die abgelegte ZIP-Datei beantwortet haben, denke ich, dass sie die folgende Datenstruktur hat.
Entpacktes Verzeichnis |-- FG-GML-0000-00-DEM10B.zip |-- FG-GML-0000-00-DEM10B.zip |-- etc...
Sie können jede ZIP-Datei weiter dekomprimieren, diesmal wird sie jedoch als ZIP-Datei behandelt, da dies etwas problematisch ist.
Der gesamte Code sieht wie folgt aus
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
#Vorbereitung
def set_value(dir, rate):
global WIDTH
global HEIGHT
#Listen und sortieren Sie Netze nach Dateinamen
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()
#Anzahl der zu verarbeitenden Daten
denominator = len(meshList)
#Anzahl der verarbeiteten Daten
numerator = 0
#Untersuchen Sie das Netz im Norden, Süden, Osten und Westen aus der Netzliste
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
#Holen Sie sich die Leinwandgröße aus dem Endnetz
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]
#Leinwandgenerierung
canvas = Image.new('RGB', (int(WIDTH / rate) * horizon, int(HEIGHT / rate) * vartical), (0, 0, 0))
#Lesen und Abbilden von Dateien
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)
#Bildgebungsmethode
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)
#Maschennummer
mesh = dem.find('xml:mesh', namespace).text
#Anzahl der angeordneten Zellen(Der tatsächliche Wert wird um 1 erhöht)
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
#Bildgrößeneinstellung(Die Anzahl der Daten==Anzahl der Pixel)
dataSize = highX * highY
#Array von Höhendaten
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])
#Datenstartkoordinaten
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
#Wenn die Anzahl der Daten nicht ausreicht(Wenn über und unter dem Bild Leerzeichen vorhanden sind)
if(len(dataNp) < dataSize):
add = []
#Wenn die Daten am unteren Rand des Bildes nicht ausreichen
if(startPosition == 0):
for i in range(dataSize - len(dataNp)):
add.append(0)
dataNp.extend(add)
#Wenn die Daten oben und unten im Bild nicht ausreichen
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)
#8-Bit-Ganzzahlkonvertierung von Daten
dataNp = (dataNp / 15).astype(np.uint8) #Teilen Sie dataNp durch 15, um den höchsten Punkt des Berges Fuji bei 255 zu erreichen
data = dataNp.reshape(highY, highX)
#Abbildung numerischer Höhendaten
pilImg = Image.fromarray(np.uint8(data))
pilImg = pilImg.resize((int(highX / rate), int(highY / rate)), Image.LANCZOS) # NEAREST
#Berechnen Sie die Koordinaten des einzufügenden Netzes
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)))
#Geben Sie das Verzeichnis ein, das die ZIP-Datei enthält
DIR = sys.argv[1]
RATE = int(sys.argv[2])
set_value(DIR, RATE)
Wenn Sie Python ausführen, sollte es funktionieren, wenn Sie ein Verzeichnis mit ZIP-Dateien angeben
Ergänzung
Wo wird eine Zeichenfläche mit der Methode set_value generiert?
canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))
Cocos 1125 und 750 schreiben jeweils nur manuell die Anzahl der Zellen in die numerischen Höhendaten.
Da die Werte selbst mit den im Code angezeigten Werten für highX und highY identisch sind, ist dies zweimal problematisch, und highX und Y suchen bereits jedes Mal nach Werten, wenn ein Bild generiert wird, und es handelt sich bereits um viel Scheißcode.
(Ich werde es vielleicht noch einmal umschreiben)
Ich bin froh, wenn jemand es umgestaltet ...
Und
dataNp = (dataNp / 15).astype(np.uint8)
Was ist 15? Die Antwort finden Sie im vorherigen Artikel.
https://qiita.com/artistan/items/99266407702d4152e9d5
Ändern wir die Größe des großen Bildes
Wie ich am Anfang geschrieben habe, ist es unpraktisch, mit dem Code so wie er ist ein großes Bild zu erzeugen.
Also der folgende Teil der set_value-Methode
canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))
Von
canvas = Image.new('RGB', (1125/10 * horizon, 750/10 * vartical), (0, 0, 0))
Und ändern Sie die Größe des auskommentierten pilImage
pilImg = pilImg.resize((int(highX)/ 10, int(highY) / 10), Image.LANCZOS) # NEAREST
Wenn Sie dies tun, ist es praktisch, ein auf 1/10 reduziertes Bild zu erzeugen.
Endlich
Hier ist die Präfektur Osaka, die mit dem obigen Code generiert und entsprechend reduziert wurde.
Osaka ist nicht interessant, weil es nicht so hoch ist, aber ich denke, es ist sehenswert, ob es ein rauerer Ort ist.
Auf einer Grauskala sieht die Form des Berges aus wie eine Blattader, die irgendwie frisch ist.
Ich konnte es nur unverständlich erklären, aber danke fürs Lesen.
Recommended Posts