In [1]:
import json
In [2]:
# Load cat data
db = json.load(open('cat-position-data.json'))
In [15]:
# Print scale factors of map images to get to 1 pixel per metre
# This is used to resize each layer in road.xcf which is 1PPM.
for cat in db['cats']:
    print('{0}: {1:.2f}%'.format(cat['name'], 100.0/cat['pixelsPerMeter']))
Ginger: 67.96%
Chip: 33.98%
Sooty: 22.65%
Orlando: 67.96%
Hermie: 47.97%
Phoebe: 33.98%
Deebee: 33.98%
Kato: 33.98%
Coco: 50.97%
Rosie: 25.48%
In [17]:
from PIL import Image
In [78]:
# Get the size of each map image in pixels and hence pixels in road.tiff
image_sizes = {}
image_scales = {}
for cat in db['cats']:
    im = Image.open(cat['name'].lower() + '.jpg')
    w, h = im.size
    scale = 1.0/cat['pixelsPerMeter']
    w = int(w*scale); h = int(h * scale)
    print('{0}: {1}x{2}'.format(cat['name'], w, h))
    image_sizes[cat['name'].lower()] = (w,h)
    image_scales[cat['name'].lower()] = scale
Ginger: 407x407
Chip: 203x203
Sooty: 135x135
Orlando: 407x407
Hermie: 287x287
Phoebe: 203x203
Deebee: 203x203
Kato: 203x203
Coco: 305x305
Rosie: 152x152
In [31]:
# This is a list of the X, Y positions of the top-left (i.e. minimum) pixel co-ord of each image in road.xcf
image_locs = {
  'chip': (2159, 2034),
  'coco': (1636, 2007),
  'deebee': (1807, 2038),
  'ginger': (1294, 2174),
  'hermie': (1619, 1880),
  'kato': (1948, 2348),
  'phoebe': (1948, 2348),
  'orlando': (0, 0), # not geo-referenced
  'rosie': (1870, 2165),
  'sooty': (1481, 2051),
}
In [33]:
from osgeo import gdal
In [39]:
# Open road map and get the geo transform
map_ds = gdal.Open('road.tiff')
ox, xs, _, oy, _, ys = map_ds.GetGeoTransform()
def pixel_to_grid(x, y):
    return (ox + xs * x, oy + ys * y)
In [40]:
pixel_to_grid(0,0)
Out[40]:
(501000.0, 146000.0)
In [41]:
pixel_to_grid(2159, 2034)
Out[41]:
(503159.0, 143966.0)
In [54]:
from osgeo import osr
In [72]:
# Setup an OS grid -> WGS 84 projection
bng = osr.SpatialReference()
bng.ImportFromEPSG(27700)
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
Out[72]:
0
In [73]:
bng2wgs84 = osr.CoordinateTransformation(bng, wgs84)
In [86]:
# Project the points for each cat into os grid and thence WGS84
for cat in db['cats']:
    name = cat['name'].lower()
    dx, dy = image_locs[name]
    
    # sckip unlocated images
    if dx == 0 and dy == 0:
        continue
    
    scale = image_scales[name]
    
    coords = []
    for pt in cat['points']:
        if not 'projected' in pt:
            continue
        x, y = tuple(float(v) for v in pt['projected'].split(','))
        x *= scale; y *= scale
        gx, gy = pixel_to_grid(x+dx, y+dy)
        
        # To WGS84
        x, y, _ = bng2wgs84.TransformPoint(gx, gy)
        coords.append((x, y))

    geometry = { 'type': 'LineString', 'coordinates': coords }
    feature = { 'type': 'Feature', 'properties': cat, 'geometry': geometry }
    json.dump(feature, open(name + '.geojson', 'w'))
In [ ]: