Try to image the elevation data of the Geographical Survey Institute with Python
A program to convert the digital elevation model of the Geographical Survey Institute to PNG using Python created by my superior's hobby () during my previous job.
It's a waste to put it to sleep, so I'll leave it to the next learner
Needless to say, the code below is not the "best code"!
Also, I'm not a Python expert, so I'm writing while thinking, "I'd be happy if someone refactored it."
* Supplement (2020/02/14) This code is bad because it was written by me in the first year. Please be careful.
What you want to do
I want to get the numerical elevation data of the Geographical Survey Institute and make it a grayscale image using Python
We will generate one PNG file for each data file
Please refer to the following Qiita article for the explanation of the xml file that can be dropped because it is very easy to understand.
https://qiita.com/tobira-code/items/43a23362f356198adce2
Or rather, with this article, I don't have to write a new article ...
Prepare the data
Immediately, select any mesh from the download page of the Geographical Survey Institute and get it.
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.
Membership registration is required to download the data!
This time, we will use the data of "Topographic map contour lines" of "10m mesh", so select the radio button and check button on the left side as such, and click any mesh.
From the dropped Zip file, unzip the file "FG-GML-0000-0-dem10b-yyyymmdd.xml" and place it in an appropriate directory.
This completes the data preparation.
The code is below.
import os
import re
import sys
import numpy as np
from PIL import Image
import xml.etree.ElementTree as ET
#xml file input
DATA = sys.argv[1]
tree = ET.parse(DATA)
root = tree.getroot()
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)
imgSize = 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(imgSize)
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
##Start image generation
#When the number of data is insufficient(If there is space above and below the image)
if(len(dataNp) < imgSize):
add = []
#When the data at the bottom of the image is insufficient
if(startPosition == 0):
for i in range(imgSize - 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(imgSize - 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), int(highY)), Image.LANCZOS) # NEAREST
canvas = Image.new('RGB', (highX, highY), (0, 0, 0))
canvas.paste(pilImg, (0, 0))
canvas.save('dem.png', 'PNG', quality=100)
Execution is OK with python .py file .xml file .
Supplement
Search for various data
You can use the above method to get the first mesh number, cell array number, elevation data, etc. You can look for it with a regular expression by turning the loop like .
For me, I prefer to search by turning the loop, so I recommend it because it's easier, but then I'll cover this article ...
"DataNp / 15" that appears in the second half
I'm sorry for the gorigori magic number, but it's kind of annoying, so I leave it as it is.
In the digital elevation model, the elevation of each point except the water surface is stored as it is as a numerical value (although I do not know if it is the case for all DEMs). The point at the top of Mt. Fuji should be stored in an xml file such as 3776.0.
If that value is directly converted to a grayscale PNG image, the black and white image will be represented by a value from 0 to 255, so all lands with an altitude of 255 m or higher will be completely white. In order to prevent this, the above code is adjusted by dividing by 15 in order to make it 255 or less based on the altitude of 3770m at the top of Mt. Fuji.
So, if you want to express Mt. Hakodate, which has an altitude of 333m, from 0 to 255 instead of the summit of Mt. Fuji, you should divide it by 1.3.
Here, it seems easy to check the maximum value of the xml file and set the value to be divided automatically.
Finally
Like this, when I converted Hakodate to PNG, it became like this.
The altitude of the city is almost 0, and I often feel that the altitude of Mt. Hakodate suddenly rises.
However, isn't it awkward?
So, next time, I will write an article that reads multiple files, sorts them by mesh number, and makes a big picture.
This will give you a pretty interesting picture.
I think this article should be reviewed in the near future and rewritten to make it easier to understand. .. .. (Maybe not)
Recommended Posts