Introduction

Describes how to use the digital elevation model (DEM) published by the Geospatial Information Authority of Japan. Describes how to convert DEM data to 2D elevation data.

file name

Base map information (digital elevation model) data 5m mesh (elevation) The file name of is in the following format.

format

FG-GML-pppp-qq-rr-DEM5X-yyyymmdd.xml

pppp is the primary mesh number qq is the secondary mesh number rr is the tertiary mesh number X is a surveying method, A: laser surveying, B: photogrammetry. yyyymmdd is the date.

--Example

FG-GML-5338-06-89-DEM5A-20130702.xml

Mesh number

The primary mesh is a rectangle with a width of 1 degree longitude and a height of 40 minutes latitude. The primary mesh number is represented by the latitude and longitude of the southwestern end of this rectangle. The primary mesh number is represented by 4 digits. The first two digits are 1.5 times the latitude. The last two digits are the last two digits of longitude.

The secondary mesh is the primary mesh divided vertically and horizontally into eight equal parts. Number (north-south, east-west) with the southwestern end as (0,0). The range of numbers is 0-7. The secondary mesh number is a two-digit number arranged in the order of north-south and east-west.

The tertiary mesh is the secondary mesh divided vertically and horizontally into 10 equal parts. Numbering is the same idea as the secondary mesh.

File contents

Only the main information is extracted

--mesh number

It is an array of primary, secondary, and tertiary mesh numbers.

  <mesh>53380689</mesh>

GridEnvelope

Describes the number of grid cell arrays that make up the mesh. The low attribute is a number indicating the northwestern end. It is always (0,0). The high attribute is a number indicating the southeastern end. The 5m mesh is fixed at (224,149) and the 10m mesh is fixed at (1124,749).

  <gml:limits>
  	<gml:GridEnvelope>
  		<gml:low>0 0</gml:low>
  		<gml:high>224 149</gml:high>
  	</gml:GridEnvelope>
  </gml:limits>

DataBlock

Elevation data for each grid cell is listed. The data enclosed in the tupleList tag is the elevation data. If there is no water surface or elevation data, the elevation value will be "-9999.".

  <gml:DataBlock>
  	<gml:rangeParameters>
<gml: QuantityList uom = "DEM configuration point" /> </gml:rangeParameters> gml:tupleList (Elevation data) Example) Ground surface, 1055.33 </gml:tupleList> </gml:DataBlock>

coverageFunction

  			<gml:sequenceRule  order="+x-y">Linear</gml:sequenceRule>

Specify the sequence of data with sequenceRule. The "+ x" in order indicates that x is in the positive direction (west → east), and the "-y" indicates that y is in the negative direction (north → south). From the northwest end to the southeast end, the data arranged in a row as shown below is placed in the tupleList tag.

  <gml:startPoint>173 0</gml:startPoint>

Specify the start point of the grid cell with startPoint. ** The data placed in the DataBlock will be the grid cell numbers (173,0) to (224,149). ** ** ** Please note that the actual number of data is not necessarily (gml: low + 1) * (gml: high + 1). ** **

Data retrieval

Now that we know the contents of the XML, we will extract the DEM data as 2D elevation data. I'll use Python to retrieve it. It is a 2D ndarray type of Numpy. The array size is fixed to the number of columns (gml: low + 1) and the number of rows (gml: high + 1). All data before startPoint is set to "-9999.".

`hyoukou.py`


import sys
import re
import numpy as np
import matplotlib.pyplot as plt

with open(sys.argv[1], "r") as f:
    # mesh
    r = re.compile("<mesh>(.+)</mesh>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            mesh = m.group(1)
            break
    # grid len
    r = re.compile("<gml:high>(.+) (.+)</gml:high>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            xlen = int(m.group(1)) + 1
            ylen = int(m.group(2)) + 1
            break
    # start
    r = re.compile("<gml:tupleList>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            break
    # data
    data = []
    r  = re.compile("</gml:tupleList>")
    r2 = re.compile(",(.+)")
    for ln in f:
        m = r.search(ln)
        if m != None:
            break
        else:
            m = r2.search(ln)
            if m != None:
                data.append(float(m.group(1)))
    # start point
    startx = starty = 0
    r = re.compile("<gml:startPoint>(.+) (.+)</gml:startPoint>")
    for ln in f:
        m = r.search(ln)
        if m != None:
            startx = int(m.group(1))
            starty = int(m.group(2))
            break

data2 = np.empty(xlen*ylen)
start_pos = starty*xlen + startx

for i in range(xlen*ylen):
    if i < start_pos:
        data2[i] = -9999.
    else:
        data2[i] = data[i-start_pos]

data = data2.reshape(ylen, xlen)

#print(data)
plt.imshow(data)
plt.colorbar()
plt.show()

How to call hyoukou.py Call with the DEM file name as an argument. Display 2D elevation data as an image using matplotlib. If there is no data loss (-9999.), It can be displayed normally.

>python hyoukou.py FG-GML-5338-06-89-DEM5B-20130702.xml

[PYTHON] Basic map information-Use digital elevation model