import json
# Load cat data
db = json.load(open('cat-position-data.json'))
# 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%
from PIL import Image
# 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
# 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),
}
from osgeo import gdal
# 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)
pixel_to_grid(0,0)
(501000.0, 146000.0)
pixel_to_grid(2159, 2034)
(503159.0, 143966.0)
from osgeo import osr
# Setup an OS grid -> WGS 84 projection
bng = osr.SpatialReference()
bng.ImportFromEPSG(27700)
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
0
bng2wgs84 = osr.CoordinateTransformation(bng, wgs84)
# 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'))