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.

ichimai

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

Try to image the elevation data of the Geographical Survey Institute with Python
I tried to find the entropy of the image with python
Try scraping the data of COVID-19 in Tokyo with Python
Try to automate the operation of network devices with Python
Try to extract the features of the sensor data with CNN
[Python] Try to graph from the image of Ring Fit [OCR]
Try to solve the shortest path with Python + NetworkX + social data
Return the image data with Flask of Python and draw it to the canvas element of HTML
How to crop the lower right part of the image with Python OpenCV
Get images of OpenStreetMap and Geographical Survey Institute maps with Python + py-staticmaps
Try to solve the man-machine chart with Python
Get images of OpenStreetMap and Geographical Survey Institute maps with Python + staticmap
[Introduction to Python] How to get the index of data with a for statement
Try to solve the programming challenge book with python3
How to scrape image data from flickr with python
Convert the image in .zip to PDF with Python
Try to get the contents of Word with Golang
Extract the band information of raster data with python
Try to measure the position of the object on the desk (real coordinate system) from the camera image with Python + OpenCV
Try to create a battle record table with matplotlib from the data of "Schedule-kun"
Try to get the function list of Python> os package
The story of rubyist struggling with python :: Dict data with pycall
[Homology] Count the number of holes in data with Python
Try to decipher the garbled attachment file name with Python
[Introduction to Python] How to get data with the listdir function
Get the source of the page to load infinitely with python.
Try to operate Facebook with Python
Try blurring the image with opencv2
Try to import to the database by manipulating ShapeFile of national land numerical information with Python
Try to visualize the nutrients of corn flakes that M-1 champion Milkboy said with Python
Save the results of crawling with Scrapy to the Google Data Store
[Cloudian # 9] Try to display the metadata of the object in Python (boto3)
Practical exercise of data analysis with Python ~ 2016 New Coder Survey Edition ~
First python ② Try to write code while examining the features of python
Try to solve the N Queens problem with SA of PyQUBO
I want to output the beginning of the next month with Python
Output the contents of ~ .xlsx in the folder to HTML with Python
Try to extract a character string from an image with Python3
Consider the speed of processing to shift the image buffer with numpy.ndarray
From the introduction of JUMAN ++ to morphological analysis of Japanese with Python
I tried to improve the efficiency of daily work with Python
Try to get CloudWatch metrics with re: dash python data source
PhytoMine-I tried to get the genetic information of plants with Python
Try to calculate the position of the transmitter from the radio wave propagation model with python [Wi-Fi, Beacon]
Basics of binarized image processing with Python
Try logging in to qiita with Python
HTML email with image to send with python
Try working with binary data in Python
Check the existence of the file with python
Introduction to Python Image Inflating Image inflating with ImageDataGenerator
Convert Excel data to JSON with python
Convert FX 1-minute data to 5-minute data with Python
[Python] Try pydash of the Python version of lodash
Try converting to tidy data with pandas
Python amateurs try to summarize the list ①
Drawing with Matrix-Reinventor of Python Image Processing-
Recommendation of Altair! Data visualization with Python
Try to generate an image with aliasing
[Introduction to Data Scientists] Basics of Python ♬
The road to compiling to Python 3 with Thrift
Sample to convert image to Wavelet with Python