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

Working with Digital Elevation Models (DEM) is an interesting extension to creating maps, which are usually a flat representation of part of the world. I really want to find a way to show elevation in a way that is a bit different from a flat map. Working with the OS DEM data has whet my appetite to try something new, but first I need a map to work with.

The most detailed DEM data I have is based on Ordnance Survey OpenData, so creating a map in the OS projection will be useful. I use TileMill to create maps from OSM data.

Firstly I needed OSM data in the Ordnance Survey projection. That means loading some OSM data into a fresh PostgreSQL database. I created a PostgreSQL database and, as usual, add the extension for PostGIS. this creates a table called spacial_ref_sys that includes the OS projection, amongst many others. I often add the hstore extension too, but this is a simple map so I didn't need it.

createdb -E UTF8 EYOS
echo "CREATE EXTENSION postgis;" | psql -d EYOS
I loaded an extract of OSM data using the usual osm2pgsql utility except the projection was needed too to convert the data to OS projection as it is loaded.

osm2pgsql --slim -d EYOS -C 1024 ey.osm.pbf --proj 27700
I decided to add a coastline, so that needed to be in OS projection too. OSM coastlines are handled differently from all other data. They are extracted from the main DB, checked for consistency and created into a shapefile for the world. This is known as processed_p.shp. I have my own copy with a cut-down version with only the British Isles in it to make rendering a bit quicker. I reprojected that to a copy in OS projection using OGR2OGR, part of the Geospatial Data Abstraction Library



ogr2ogr -t_srs 'EPSG:27700' -s_srs 'EPSG:3857' coast_bi_os.shp coast_bi.shp

Armed with all of this I could now start TileMill and add the layers I need for the map. Each of the layers, including the coast shapefile, needed a custom projection. This is:
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy  +units=m +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +nodefs
 I got this from the PostgreSQL table postgis created above. Once I had designed the map as I wanted it I exported the Mapnik XML and ran it through Mapnik. I discovered that the Mapnik XML was not quite right. It needed to have the third line changed so the srs part matches the custom projection above. There doesn't seem to be a way to set this in TileMill, so a manual edit was needed.

To run the Mapnik XML through Mapnik I used the following python code:

#!/usr/bin/python

# generate a map image in OS projection epsg:27700

import mapnik
import sys, os
def drawMap(filename, west,south,east,north):
    print(filename)
    sz = 5000
    ossrs = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"
    m = mapnik.Map(sz,sz,ossrs)
    mapnik.load_map(m,"osproj.xml")
    bbox = mapnik.Envelope(west,south,east,north)
    m.zoom_to_box(bbox)
    im = mapnik.Image(sz,sz)
    mapnik.render(m, im)
    view = im.view(0,0,sz,sz) # x,y,width,height
    view.save(filename,'png')

if __name__ == '__main__':
    drawMap('cott.png',500000,430000,510000,440000)
 

 That long-winded projection was needed again. Notice the coordinates in the drawMap function are OS coordinates, not longitude and latitude. Everything must match the chosen projection.

This gives me an image of the map in the OS projection, but the style could be any style you choose, though I'd be wary of copying the OS style too closely. This will now match the DEM data if they are combined. My style is still a bit stark and only renders a few objects, but it is something to work with.

Next I need to use it imaginatively.

Friday 8 May 2015

Heights

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.

SRTM
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__":
    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 pt                    lowest=pt
                if pt>highest:
                    highest=pt
                im.putpixel((e,n),getpix(pt))
       
    print 'lowest:{0}, highest:{1}'.format(lowest,highest)
    im.save('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):
    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')