Thursday, 21 May 2015

Not Alone

I went out looking for new developments in East Yorkshire today. Looks like I was not alone.

Monday, 11 May 2015

Heights & OS




createdb -E UTF8 EYOS
echo "CREATE EXTENSION postgis;" | psql -d EYOS






# generate a map image in OS projection epsg:27700

import mapnik
import sys, os
    sz = 5000
    m = mapnik.Map(sz,sz,ossrs)
    bbox = mapnik.Envelope(west,south,east,north)
    im = mapnik.Image(sz,sz)
    mapnik.render(m, im)

if __name__ == '__main__':



Next I need to use it imaginatively.

Friday, 8 May 2015


I've been working on something locally for a while that benefits from maps. It needs height information displayed so I thought I'd take a closer look at what was available, especially Digital Elevation Model (DEM) data

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.

The SRTM is published as a 1° square. I read the height values and displayed them as a shade of green since human eyes can distinguish more share variations in green than any other colour. Any voids I show as black There were nine pixels in this square) and any value with a small negative value (small so not a void) I show as blue. There's a lot of interest in this which I'd not noticed looking at contours. The dark area top left is the Vale of York, the green area top centre is the bottom end of the Yorkshire Wolds You can just make out the Humber estuary just above the centre and to the right. The bright green area bottom right is part of the Lincolnshire Wolds. The valleys with tributaries feeding into the Vale of York are interesting. None of those exist as rivers or streams today, so I expect they are remnants of the retreating ice caps about ten thousand years ago when the ground was still permafrost so any melting water cut river channels. Today the water table is much lower with the chalk of the Wolds allowing water to drain into it.

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
I created a similar, image from the OS DEM data. I deliberately emphasised the height differences more than the SRTM image. The area is much smaller than the SRTM area but more detailed. bright green on the left is the edge of the Yorkshire Wolds. The blue line is the river Hull which cuts through the middle of the city of Hull. For comparison the bottom of the the blue smudge on the SRTM image is approximately where the OS image is. Again valleys are shown, though this time running west to east. Again they are dry (though very occasionally not which is part of what I'm investigating). I've decided that there is enough detail in the OS area and that it is big enough, perhaps with one more alongside it, to show what I want. so I'll work with that.

More of what to do with it later.

The python code to produce the SRTM image is here:
#!/usr/bin/env python
# -*- 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__":
    tile = "N53W001.hgt"
    #make the new empty (white) image
    im ="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 =
                hite=struct.unpack('>h', buf)
                #print '{0} {1} {2}'.format(n,e,hite[0])
                if pt == -32768:
                    print 'VOID {0} {1} {2}'.format(n,e,pt)
                if pt                    lowest=pt
                if pt>highest:
    print 'lowest:{0}, highest:{1}'.format(lowest,highest)'h.png')
The code to produce the OS image is here:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import struct
from PIL import Image

def getpix(num):
    if val < 0:
        return (0,0,255)
    return (0,int(val)*3,0)

if __name__ == "__main__":
    osf = "TA03.asc"
    #make the new empty (white) image
    im ="RGB", (200, 200), "white")
    with open(osf, "rb") as f:
        for i in range(5,205):
            s=lines[i].split(' ')
            for idx, val in enumerate(s):