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.


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.

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
|-- 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

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])

    #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:
        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 ='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))'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 as zfile:
            root = ET.fromstring(
            namespace = {
                'xml': '',
                'gml': ''
            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
                    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)):
                #When the data at the top and bottom of the image is insufficient
                    for i in range(startPosition):
                    dataNp[0:0] = add
                    add = []
                    for i in range(dataSize - len(dataNp) - len(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


Where to generate canvas with set_value method

canvas ='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 ...

dataNp = (dataNp / 15).astype(np.uint8)

What's 15? Please refer to the previous article for the answer.

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 ='RGB', (1125 * horizon, 750 * vartical), (0, 0, 0))


canvas ='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.


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

Make multiple numerical elevation data into a single picture with Python
Make holiday data into a data frame with pandas
Make a fortune with Python
Let's make a GUI with python.
Make a recommender system with python
Let's make a graph with python! !!
Add a Python data source with Redash
Let's make a voice slowly with Python
Let's make a web framework with Python! (1)
Make a desktop app with Python with Electron
Let's make a Twitter Bot with Python!
Let's make a web framework with Python! (2)
Make a decision tree from 0 with Python and understand it (4. Data structure)
Make Python scripts into Windows-executable .exes with Pyinstaller
Make a Twitter trend bot with heroku + Python
[Python] Make a game with Pyxel-Use an editor-
I want to make a game with Python
Try to make a "cryptanalysis" cipher with Python
A story stuck with handling Python binary data
[Python] Make a simple maze game with Pyxel
Folium: Visualize data on a map with Python
Let's replace UWSC with Python (5) Let's make a Robot
Try to make a dihedral group with Python
Make JSON into CSV with Python from Splunk
Extract data from a web page with Python
Make one repeating string with a Python regular expression.
Data analysis with python 2
Get financial data with python (then a little tinkering)
[Practice] Make a Watson app with Python! # 2 [Translation function]
[Practice] Make a Watson app with Python! # 1 [Language discrimination]
Make a simple Slackbot with interactive button in python
[Let's play with Python] Make a household account book
Let's make a simple game with Python 3 and iPhone
Make a breakpoint on the c layer with python
Fill the background with a single color with OpenCV2 + Python
Until you insert data into a spreadsheet in Python
A server that echoes data POSTed with flask / python
Make a CSV formatting tool with Python Pandas PyInstaller
What is God? Make a simple chatbot with python
[Super easy] Let's make a LINE BOT with Python.
Data analysis with Python
A Python program that converts ical data into text
Numerical calculation with Python
Generate multiple HTML files at once by pouring JSON data into an HTML template with Python
Let's make a websocket client with Python. (Access token authentication)
[Practice] Make a Watson app with Python! # 3 [Natural language classification]
Embed a Python interpreter into a C ++ app with pybind11 + cmake
I tried to make various "dummy data" with Python faker
Associate Python Enum with a function and make it Callable
Experiment to make a self-catering PDF for Kindle with Python
Create applications, register data, and share with a single email
Stylish technique for pasting CSV data into Excel with Python
[Python] Make AWS resources mocked with Moto into pytest fixtures
[Python] Make a simple maze game with Pyxel-Make enemies appear-
Sample data created with python
Get Youtube data with python
Make Puyo Puyo AI with Python
Make a bookmarklet in Python
Create a directory with python
Make a fire with kdeplot
Read json data with python