Informations sur la carte de base à l'aide de la conversion Geotiff Python des données numériques d'élévation
introduction
Les données du titre sont diffusées sur ce site. Des données d'altitude (DEM) très détaillées et précises au Japon sont disponibles. Ce qui est téléchargé est un fichier au format zip pour chaque maillage de zone de division secondaire défini par le Bureau des statistiques du Ministère de l'intérieur et des communications. Il contient des données d'élévation au format xml. Le maillage régional sera décrit plus loin.
Le code se trouve à la fin de cet article. Un exemple est le suivant.
** (Étant donné que le nom du fichier est également utilisé comme méta-information, veuillez l'appliquer à celui téléchargé) **
commander
python zdem2tif.py FG-GML-1234-56-DEM10B.zip FG-GML-6543-21-DEM5A.zip ...
Une image Geotiff est sortie pour chaque fichier zip (= maillage de la zone de partition secondaire). Prend en charge les résolutions de 5 m et 10 m. Nous avons confirmé l'opération avec la série python3. Le package gdal doit être installé.
Cette fois, j'ai beaucoup fait référence au code placé dans ce site. Cependant, s'il est laissé tel quel, un bogue se produira en cas de perte de données, donc ce code a été corrigé.
Il semble y avoir divers autres outils, mais j'ai trouvé ce code car il peut être exécuté sur des commandes, il peut être partagé avec des étrangers car il est écrit en anglais et il n'est pas nécessaire d'étendre le zip. ..
De plus, comme je l'ai remarqué après la publication, qiita a déjà un article sur les informations de base sur la carte et python. L'explication du format de fichier est plus détaillée ici.
À propos du maillage régional
La définition détaillée est expliquée dans ce document. Les maillages régionaux sont définis dans la 1ère à la 3ème division en fonction de la résolution spatiale. Le maillage principal reçoit un code de partition de 4 chiffres, le maillage secondaire reçoit 2 chiffres supplémentaires et le maillage tertiaire se voit attribuer 2 chiffres supplémentaires. La latitude et la longitude des quatre coins du compartiment peuvent être connues à l'aide du numéro du code de compartiment.
Description du code source
À l'avenir, il peut y avoir des changements dans les spécifications des données distribuées, et il est possible qu'elles ne soient pas prises en charge dans les zones que l'auteur n'a pas confirmées, je vais donc expliquer brièvement le code.
L'ensemble du flux (convertir) est
- Obtenez des informations de maillage de zone et la résolution DEM à partir du nom de fichier zip
--Définissez un tableau de fichiers de sortie (jusqu'ici
DEMOutArray,mesh_info_from_zipname) - Boucle sur le fichier xml dans le fichier zip (
pour le nom de fichier dans la liste de fichiers) - Récupérez les informations requises à partir du fichier xml (
XMLContents,read_xml) --Générer un tableau DEM pour chaque fichier xml (DEMInArray) - Écraser le tableau de sortie (ʻupdate_raster
) --Enregistrer au format Tiff (save_raster`)
Obtenez des informations de maillage de zone et la résolution DEM à partir du nom de fichier zip
Calculez la latitude et la longitude des quatre coins en lisant le code de partition principale à 4 chiffres et le code de partition secondaire à 2 chiffres. En outre, la résolution est interprétée selon que le nom du fichier est écrit DEM5X ou DEM10X. Pour chaque résolution, la taille du tableau d'un fichier xml est fixe (par exemple, dans le cas de 10mDEM, on suppose qu'il n'y a qu'un seul xml dans le zip).
Définir un tableau de fichiers de sortie
Créez à l'avance la taille d'un compartiment secondaire. Dans le cas de 5mDEM, on suppose qu'un maximum de 10 x 10 xml (partition tertiaire) est contenu dans un zip (partition secondaire). Cependant, il semble qu'il y ait de nombreuses zones qui n'ont pas été observées, et il y a des cas où seulement environ 10 xml existent. Une valeur aberrante de -9999. A été utilisée pour l'emplacement correspondant à la zone où xml n'existe pas.
Boucle sur le fichier xml dans un fichier zip
J'ai d'abord découvert le module qui est inclus par défaut appelé zipfile.
En utilisant cela, vous pouvez lire les fichiers dans le zip sans décompresser chacun d'eux.
Obtenez les informations requises à partir du fichier xml
Pour 5mDEM, lisez les deux derniers chiffres du code de partition tertiaire dans le nom de fichier pour déterminer la position d'insertion dans le tableau de sortie. Dans le cas de 10mDEM, comme mentionné ci-dessus, on suppose qu'il n'y a pas de description du code de partition tertiaire et qu'il n'y a qu'un seul xml. Obtenez les informations nécessaires au format texte pour le moment. Les règles pour les balises en xml, etc. sont supposées comme dans le code. Si vous lisez xml avec des règles différentes, cela devrait apparaître immédiatement comme un bogue.
Générer un tableau DEM pour chaque fichier xml
Générez un tableau en utilisant les informations décrites dans le fichier xml. Spécifications du fichier a une définition relative aux données DEM. Selon lui, la division terrestre pour chaque pixel est soit {surface du sol, surface de surface, surface de l'eau de mer, surface de l'eau intérieure, pas de données, etc.}, et le nombre de -9999 est entré pour les endroits où la valeur n'est pas définie. ..
Dans le code référencé à l'origine, les informations DEM sont acquises selon qu'il s'agit de "surface au sol" ou non, mais dans ce code, elles sont divisées selon que la valeur DEM en xml est -9999 ou non (indiquée par "Autre"). Parce qu'il y avait une valeur DEM dans la zone. La raison de dire "autre" est inconnue). Pour les endroits où DEM n'est pas défini, la valeur d'écart -1 est appliquée, et le cas où le fichier xml mentionné ci-dessus n'est pas en premier lieu et l'utilisation de la valeur d'écart sont séparés.
Écraser le tableau de sortie
En ce qui concerne l'orientation du réseau, on suppose que la direction de la latitude est définie comme vers le nord. Il y a une donnée obtenue à partir de xml appelée gml: sequenceRule, et le fait qu'elle soit écrite comme '+ x-y' le supporte. Si c'est `` + x + y '' etc., c'est probablement face au sud, mais je n'ai pas divisé les conditions autour de cela.
Enregistrer au format Tiff
Il n'y a rien de spécial à dire.
Dessin du tiff de sortie
Fichiers convertis téléchargés avec 5 mDEM en haut et 10 mDEM en bas pour la même zone de tracé secondaire (valeurs d'axe presque inédites, etc.). Les valeurs négatives sont masquées. On constate que le taux d'abondance des données à chaque résolution diffère selon la région. Je serais heureux si je le savais avant de le télécharger, mais ne manque-t-il aucune information? .. ..
Code source complet
Enfin tout le 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 == 'Surface au sol':
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