OSM doesn't hold much height information, so when people want to display heights they turn to outside information. One such source is Shuttle Radar Topography Mission data or SRTM. One Space Shuttle mission flew around the world and mapped the heights of the ground below using radar. This data has been published as open data. It is 1 second of arc data points for the USA and 3 seconds of arc data points for the rest of the world. This gives a height data point about every 90m for the UK. There are issues with this data with some places having voids where the radar return didn't register. It is usual for people who use this with OSM to render this data as contour lines or as hill shading or both as a way of visualising the height. I thought I'd do some simple processing to be sure I understood the data format.
SRTM |
Next I looked at Ordnance Survey (OS) OpenData. They release height data as contours and spot heights in shape file format and DEM data too. The DEM is 50m spacing and should be free of voids. They use their own projection (EPSG:27700) for all of their data and this works better for the UK for some jobs. OS release some of their data in parcels based on their own grid. I am interested in a section of including Cottingham, a large village west of Hull. The OS square TA03 has Cottingham in the middle of it, so that is helpful.
OS TA03 square |
More of what to do with it later.
The python code to produce the SRTM image is here:
#!/usr/bin/env pythonThe code to produce the OS image is here:
# -*- coding: utf-8 -*-
import struct
from PIL import Image
def getpix(val):
if val == -32768:
return (0,0,0)
if val < 0:
return (0,0,255)
return (0,int(val),0)
if __name__ == "__main__":
top=54
left=5
highest=-5000
lowest=5000
tile = "N53W001.hgt"
#make the new empty (white) image
im = Image.new("RGB", (1201, 1201), "white")
with open(tile, "rb") as f:
#print get_sample(tile, n, e)
# scan through each of the heights in the file and colour a pixel
for n in range(1201):
for e in range(1201):
buf = f.read(2)
hite=struct.unpack('>h', buf)
#print '{0} {1} {2}'.format(n,e,hite[0])
pt=hite[0]
if pt == -32768:
print 'VOID {0} {1} {2}'.format(n,e,pt)
if ptlowest=pt
if pt>highest:
highest=pt
im.putpixel((e,n),getpix(pt))
print 'lowest:{0}, highest:{1}'.format(lowest,highest)
im.save('h.png')
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import struct
from PIL import Image
def getpix(num):
val=round(num)
if val < 0:
return (0,0,255)
return (0,int(val)*3,0)
if __name__ == "__main__":
top=54
left=5
highest=-5000
lowest=5000
osf = "TA03.asc"
#make the new empty (white) image
im = Image.new("RGB", (200, 200), "white")
with open(osf, "rb") as f:
lines=f.readlines()
for i in range(5,205):
s=lines[i].split(' ')
for idx, val in enumerate(s):
im.putpixel((idx,i-5),getpix(float(val)))
im.save('o.png')
No comments:
Post a Comment