Make multiple numerical elevation data into a single picture with Python
Continuation of the article I wrote earlier? Something like
"Try to image the elevation data of the Geographical Survey Institute with Python"
In this article, I'd like to throw in multiple xml files and arrange them properly to make a single picture.
Preparation
This time, we will drop the entire mesh data of any prefecture.
Select and obtain any mesh from the download page of the Geographical Survey Institute.
https://fgd.gsi.go.jp/download/menu.php
Click the "Go to file selection" button in "Basic map information digital elevation model" to jump to the download screen.
This time, I will use the data of Osaka prefecture.
Please check 10 (contour line of topographic map) of 10m mesh on the left side.
In the code written in this article, the data provided by the Geographical Survey Institute is grayscaled as it is
In other words, in the case of 10m mesh data, a 1125x750 grayscale image is generated for each xml file, and they are arranged one by one to generate one large image.
In the case of Osaka prefecture, a stupid big image of 6x1125px in width and 10x750px in height will be generated, so it is better to select a small prefecture.
If you want to reduce the number of pixels in a single picture, please play with the parts described below.
By the way, once you answer the dropped ZIP file, I think that it has the following data structure.
Unzipped directory |-- FG-GML-0000-00-DEM10B.zip |-- FG-GML-0000-00-DEM10B.zip |-- etc...
You can decompress each ZIP file further, but this time it will be treated as a ZIP file because it is a little troublesome.
The whole code looks like the following
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
#Preparation
def set_value(dir, rate):
global WIDTH
global HEIGHT
#List and sort meshes by file name
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()
#Number of data to be processed
denominator = len(meshList)
#Number of processed data
numerator = 0
#Examine the mesh at the north, south, east, and west from the mesh list
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
#Get the canvas size from the end mesh
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]
#canvas generation
canvas = Image.new('RGB', (int(WIDTH / rate) * horizon, int(HEIGHT / rate) * vartical), (0, 0, 0))
#File reading and imaging
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)
#Imaging method
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)
#Mesh number
mesh = dem.find('xml:mesh', namespace).text
#Number of cell arrays(Actual value is added by 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
#Image size setting(The number of data==Number of pixels)
dataSize = highX * highY
#Array of elevation data
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])
#Data start coordinates
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
#When the number of data is insufficient(If there is space above and below the image)
if(len(dataNp) < dataSize):
add = []
#When the data at the bottom of the image is insufficient
if(startPosition == 0):
for i in range(dataSize - len(dataNp)):
add.append(0)
dataNp.extend(add)
#When the data at the top and bottom of the image is insufficient
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 integer conversion of data
dataNp = (dataNp / 15).astype(np.uint8) #If you want to fit the highest point of Mt. Fuji at 255, divide dataNp by 15.
data = dataNp.reshape(highY, highX)
#Imaging numerical elevation data
pilImg = Image.fromarray(np.uint8(data))
pilImg = pilImg.resize((int(highX / rate), int(highY / rate)), Image.LANCZOS) # NEAREST
#Calculate the coordinates of the mesh to be pasted
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)))
#Enter the directory containing the ZIP file
DIR = sys.argv[1]
RATE = int(sys.argv[2])
set_value(DIR, RATE)
When running python, it should work if you specify a directory full of ZIP files
Supplement
Where to generate canvas with set_value method
canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))
Coco's 1125 and 750 each just manually write the number of cells in the numerical elevation data.
Since the values themselves are the same as the highX and highY that appear in the code, it is troublesome twice, and the highX and Y are already various shit codes such as searching for the value every time an image is generated.
(I'll rewrite it again, maybe)
I'm glad if someone refactors it ...
And
dataNp = (dataNp / 15).astype(np.uint8)
What's 15? Please refer to the previous article for the answer.
https://qiita.com/artistan/items/99266407702d4152e9d5
Let's resize the big image
As I wrote at the beginning, it is inconvenient to generate a big image with the code as it is.
So, the following part of the set_value method
canvas = Image.new('RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))
From
canvas = Image.new('RGB', (1125/10 * horizon, 750/10 * vartical), (0, 0, 0))
And resize the commented out pilImage
pilImg = pilImg.resize((int(highX)/ 10, int(highY) / 10), Image.LANCZOS) # NEAREST
If you do, it will be convenient to generate an image reduced to 1/10.
Finally
Here is Osaka prefecture, which was generated with the above code and reduced appropriately.
Osaka isn't interesting because it's not so high in the mountains, but I think it's worth seeing if it's a more rugged place.
In grayscale, the shape of the mountain looks like a leaf vein, which is kind of fresh.
I could only explain it incomprehensible, but thank you for reading.
Recommended Posts