Basiskarteninformationen mithilfe der Python-Geotiff-Konvertierung numerischer Höhendaten
Einführung
Die Daten im Titel werden auf dieser Site verteilt. In Japan sind sehr detaillierte und genaue Höhendaten (DEM) verfügbar. Was heruntergeladen wird, ist eine Zip-Format-Datei für jedes Netz der sekundären Abteilung, das vom Statistikbüro des Ministeriums für innere Angelegenheiten und Kommunikation definiert wird. Es enthält Höhendaten im XML-Format. Das regionale Netz wird später beschrieben.
Der Code befindet sich am Ende dieses Artikels. Ein Beispiel ist wie folgt.
** (Da der Dateiname auch als Metainformation verwendet wird, wenden Sie ihn bitte auf die heruntergeladene an.) **
Befehl
python zdem2tif.py FG-GML-1234-56-DEM10B.zip FG-GML-6543-21-DEM5A.zip ...
Für jede Zip-Datei wird ein Geotiff-Bild ausgegeben (= sekundäres Partitionsbereichsnetz). Unterstützt Auflösungen von 5 m und 10 m. Wir haben die Operation mit der Python3-Serie bestätigt. Das gdal-Paket muss installiert sein.
Dieses Mal habe ich viel auf den Code in [dieser Site] verwiesen (https://tm23forest.com/contents/python-jpgis-gml-dem-geotiff). Wenn es jedoch unverändert bleibt, tritt bei Datenverlust ein Fehler auf, sodass dieser Code behoben wurde.
Es scheint, dass es verschiedene andere Tools gibt, aber ich habe mir diesen Code ausgedacht, weil er auf Befehlen ausgeführt werden kann, mit Ausländern geteilt werden kann, weil er in Englisch geschrieben ist und keine Notwendigkeit besteht, zip zu erweitern. ..
Wie ich nach dem Posten bemerkt habe, hat qiita bereits Artikel über grundlegende Karteninformationen und Python. Die Erläuterung des Dateiformats finden Sie hier ausführlicher.
Über regionale Maschen
Die detaillierte Definition wird in diesem Dokument erläutert. Regionale Netze werden in der 1. bis 3. Abteilung entsprechend der räumlichen Auflösung definiert. Dem primären Netz wird ein Partitionscode mit 4 Stellen zugewiesen, dem sekundären Netz werden zusätzliche 2 Stellen zugewiesen, und dem tertiären Netz werden zusätzliche 2 Stellen zugewiesen. Der Breiten- und Längengrad der vier Ecken des Fachs kann anhand der Nummer des Fachcodes ermittelt werden.
Beschreibung des Quellcodes
In Zukunft kann es zu Änderungen in den Spezifikationen der verteilten Daten kommen, und es besteht die Möglichkeit, dass diese in Bereichen, die der Autor nicht bestätigt hat, nicht unterstützt werden. Daher werde ich den Code kurz erläutern.
Der gesamte Fluss (convert) ist
- Abrufen von Flächennetzinformationen und DEM-Auflösung aus dem Namen der Zip-Datei
- Definieren Sie ein Array von Ausgabedateien (bisher
DEMOutArray,mesh_info_from_zipname) --Loop über XML-Datei in Zip-Datei (für Dateinamen in Dateiliste) - Holen Sie sich die erforderlichen Informationen aus der XML-Datei (
XMLContents,read_xml) - Generiere ein DEM-Array für jede XML-Datei (
DEMInArray) --Oritrite-Ausgabearray (update_raster) - Im Tiff-Format speichern (
save_raster)
Informationen zum Flächennetz und zur DEM-Auflösung erhalten Sie vom Namen der Zip-Datei
Berechnen Sie den Breiten- und Längengrad der vier Ecken, indem Sie den 4-stelligen primären Partitionscode und den 2-stelligen sekundären Partitionscode lesen. Darüber hinaus wird die Auflösung dahingehend interpretiert, ob der Dateiname als DEM5X oder DEM10X geschrieben ist. Für jede Auflösung ist die Arraygröße einer XML-Datei festgelegt (im Fall von 10 mDEM wird beispielsweise angenommen, dass sich nur eine XML in der Zip-Datei befindet).
Definieren Sie ein Array von Ausgabedateien
Erstellen Sie im Voraus die Größe eines sekundären Fachs. Im Fall von 5mDEM wird angenommen, dass maximal 10 x 10 xml (tertiäre Partition) in einer Zip (sekundäre Partition) enthalten sind. Es scheint jedoch, dass es viele Bereiche gibt, die nicht beobachtet wurden, und es gibt Fälle, in denen nur etwa 10 xml existieren. Ein Ausreißer von -9999 wurde für den Speicherort verwendet, der dem Bereich entspricht, in dem XML nicht vorhanden ist.
Schleife über XML-Datei in Zip-Datei
Ich habe zuerst etwas über das standardmäßig enthaltene Modul namens "zipfile" erfahren. Auf diese Weise können Sie die Dateien in der Zip-Datei lesen, ohne sie zu dekomprimieren.
Holen Sie sich die erforderlichen Informationen aus der XML-Datei
Lesen Sie für 5mDEM die letzten beiden Ziffern des tertiären Partitionscodes im Dateinamen, um die Einfügeposition im Ausgabearray zu bestimmen. Im Fall von 10mDEM wird, wie oben erwähnt, angenommen, dass es überhaupt keine Beschreibung des tertiären Partitionscodes gibt und es nur eine XML gibt. Holen Sie sich vorerst die notwendigen Informationen im Textformat. Die Regeln für Tags in XML usw. werden wie im Code angenommen. Wenn Sie XML mit anderen Regeln lesen, sollte dies sofort als Fehler angezeigt werden.
Generieren Sie ein DEM-Array für jede XML-Datei
Generieren Sie ein Array mit den in der XML-Datei beschriebenen Informationen. Dateispezifikationen enthält eine Definition für DEM-Daten. Demnach ist die Landteilung für jedes Pixel entweder {Bodenoberfläche, Oberflächenoberfläche, Meeresoberfläche, Binnenwasseroberfläche, keine Daten usw.}, und die Zahl von -9999 wird für Stellen eingegeben, an denen der Wert undefiniert ist. ..
In dem ursprünglich referenzierten Code werden DEM-Informationen abhängig davon erfasst, ob es sich um "Bodenoberfläche" handelt oder nicht. In diesem Code wird jedoch danach unterteilt, ob der DEM-Wert in XML -9999 beträgt oder nicht (angezeigt durch "Andere"). Weil es in der Gegend einen DEM-Wert gab. Der Grund für die Aussage "Sonstige" ist unbekannt. Für Stellen, an denen DEM nicht definiert ist, wird der Abweichungswert -1 angewendet, und der Fall, in dem die oben genannte XML-Datei nicht an erster Stelle steht, und die Verwendung des Abweichungswerts werden getrennt.
Ausgabearray überschreiben
In Bezug auf die Ausrichtung des Arrays wird angenommen, dass die Breitengradrichtung als nach Norden definiert ist. Es gibt Daten aus xml mit dem Namen "gml: sequenceRule", und die Tatsache, dass sie als "+ x-y" geschrieben sind, unterstützt dies. Wenn dies "+ x + y" usw. ist, ist es wahrscheinlich nach Süden ausgerichtet, aber ich habe die Bedingungen dafür nicht aufgeteilt.
Im Tiff-Format speichern
Es gibt nichts Besonderes zu sagen.
Zeichnung des Ausgangstiffs
Konvertierte Dateien, die mit 5 mDEM oben und 10 mDEM unten für denselben sekundären Plotbereich heruntergeladen wurden (fast unbearbeitete Achsenwerte usw.). Negative Werte werden maskiert. Es ist ersichtlich, dass die Datenhäufigkeit bei jeder Auflösung je nach Region unterschiedlich ist. Ich würde mich freuen, wenn ich es wüsste, bevor ich es heruntergeladen habe, aber fehlen keine Informationen? .. ..
Ganzer Quellcode
Endlich der ganze Code.
zdem2tif.py
import os
import sys
import zipfile
import numpy as np
import xml.etree.ElementTree as et
import osgeo.gdal
import osgeo.osr
class XMLContents(object):
def __init__(self):
self.ns = {
'ns': 'http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema',
'gml': 'http://www.opengis.net/gml/3.2',
'xsi': 'http://www.w3.org/2001/XMLSchema-instance',
'xlink': 'http://www.w3.org/1999/xlink'
}
def read_xml(self, zf, filename, dem_mode):
if dem_mode == 'DEM10':
self.posix = 0
self.posiy = 0
elif dem_mode == 'DEM5':
_split = os.path.basename(filename).split('-')
self.posix = int(_split[4][1])
self.posiy = int(_split[4][0])
with zf.open(filename, 'r') as file:
self.text = file.read().decode('utf_8')
self.root = et.fromstring(self.text)
self.dem = self.root.find('ns:DEM', self.ns)
self.coverage = self.dem.find('ns:coverage', self.ns)
self.envelope = self.coverage.find(
'gml:boundedBy//gml:Envelope', self.ns)
self.lower = self.envelope.find('gml:lowerCorner', self.ns).text
self.upper = self.envelope.find('gml:upperCorner', self.ns).text
self.grid = self.coverage.find(
'gml:gridDomain//gml:Grid//gml:limits//gml:GridEnvelope', self.ns)
self.low = self.grid.find('gml:low', self.ns).text
self.high = self.grid.find('gml:high', self.ns).text
self.tuplelist = self.coverage.find(
'gml:rangeSet//gml:DataBlock//gml:tupleList', self.ns).text
self.gridfunc = self.coverage.find(
'gml:coverageFunction//gml:GridFunction', self.ns)
self.rule = self.gridfunc.find('gml:sequenceRule', self.ns)
self.start = self.gridfunc.find('gml:startPoint', self.ns).text
if self.rule.attrib['order'] != '+x-y':
print('warning sequence order not +x-y')
if self.rule.text != 'Linear':
print('warning sequence rule not Linear')
file.close()
class DEMInArray(object):
def __init__(self, contents):
self._noValue = -1.
self.llat, self.llon = np.array(
contents.lower.split(' '), dtype=np.float64)
self.ulat, self.ulon = np.array(
contents.upper.split(' '), dtype=np.float64)
self.lowx, self.lowy = np.array(
contents.low.split(' '), dtype=np.int)
self.highx, self.highy = np.array(
contents.high.split(' '), dtype=np.int)
self.sizex, self.sizey = self.get_size()
self.x_init, self.y_init = np.array(
contents.start.split(' '), dtype=np.int)
self.datapoints = contents.tuplelist.splitlines()
self.raster = self.get_raster()
def get_size(self):
sizex = self.highx - self.lowx + 1
sizey = self.highy - self.lowy + 1
return sizex, sizey
def get_raster(self):
_x, _y = self.x_init, self.y_init
raster = np.zeros([self.sizey, self.sizex])
raster.fill(self._noValue)
for dp in self.datapoints:
if _y > self.highy:
print('DEM datapoints are unexpectedly too long')
break
s = dp.split(',')
if len(s) == 1:
continue
desc, value = s[0], np.float64(s[1])
# if desc == 'Bodenbelag':
raster[_y, _x] = value if value > -9998 else self._noValue
_x += 1
if _x > self.highx:
_x = 0
_y += 1
return raster
class DEMOutArray(object):
def __init__(self):
self._fillValue = -9999.
def initialize_raster(self):
raster = np.zeros([self.sizey, self.sizex])
raster.fill(self._fillValue)
return raster
def get_trans(self):
# +x-y
return [self.llon, self.resx, 0, self.ulat, 0, -self.resy]
def update_raster(self, tile, posiy, posix):
xmin = posix * tile.shape[1]
ymin = posiy * tile.shape[0]
xmax = xmin + tile.shape[1]
ymax = ymin + tile.shape[0]
self.raster[ymin:ymax, xmin:xmax] = tile[::-1]
def get_size(self):
sizex = (self.highx - self.lowx + 1) * self.tilex
sizey = (self.highy - self.lowy + 1) * self.tiley
return sizex, sizey
def mesh_info_from_zipname(self, zipname):
'''
Ref: Table 4 of https://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf
'''
_split = os.path.basename(zipname).split('-')
_lat1 = int(_split[2][:2])
_lon1 = int(_split[2][2:])
_lat2 = int(_split[3][0])
_lon2 = int(_split[3][1])
self.llat = _lat1 * 0.66666666666 + _lat2 * 0.08333333333
self.llon = _lon1 + 100. + _lon2 * 0.125
self.ulat = self.llat + 0.08333333333
self.ulon = self.llon + 0.125
if _split[4][:5] == 'DEM10':
self.mode = 'DEM10'
self.tilex, self.tiley = 1, 1
self.lowx, self.lowy = 0, 0
self.highx, self.highy = 1124, 749
self.sizex, self.sizey = self.get_size()
self.resx, self.resy = 1.1111111111e-04, 1.1111111111e-04
elif _split[4][:4] == 'DEM5':
self.mode = 'DEM5'
self.tilex, self.tiley = 10, 10
self.lowx, self.lowy = 0, 0
self.highx, self.highy = 224, 149
self.sizex, self.sizey = self.get_size()
self.resx, self.resy = 5.5555555555e-05, 5.5555555555e-05
else:
raise ValueError('Unexpected definition of DEM:', _split[4])
self.raster = self.initialize_raster()
def save_raster(self, out_filename):
trans = self.get_trans()
srs = osgeo.osr.SpatialReference()
srs.ImportFromEPSG(4612)
driver = osgeo.gdal.GetDriverByName('GTiff')
output = driver.Create(out_filename, self.sizex, self.sizey,
1, osgeo.gdal.GDT_Float32)
output.GetRasterBand(1).WriteArray(self.raster)
output.GetRasterBand(1).SetNoDataValue(self._fillValue)
output.SetGeoTransform(trans)
output.SetProjection(srs.ExportToWkt())
output.FlushCache()
def convert(in_filename, out_filename):
output_array = DEMOutArray()
output_array.mesh_info_from_zipname(in_filename)
dem_mode = output_array.mode
with zipfile.ZipFile(in_filename, 'r') as zf:
filelist = zf.namelist()
for filename in filelist:
print('loading', filename)
contents = XMLContents()
contents.read_xml(zf, filename, dem_mode)
posix, posiy = contents.posix, contents.posiy
dem_tile = DEMInArray(contents)
output_array.update_raster(dem_tile.raster, posiy, posix)
zf.close()
print('save', out_filename)
output_array.save_raster(out_filename)
def main(argv):
for i in range(1, len(argv)):
convert(argv[i], argv[i].replace('.zip', '.tif'))
if __name__ == '__main__':
main(sys.argv)
Recommended Posts