This website is based on the Google API, however it should be not be too hard to transfer it to another java script map viewer. Such as OpenLayers, which would remove the requirment to have the internet, to view maps. Client side scripting is done in javascript, and server side scripting in python.
The main work so far has been in the tile cutting, which is the process of making small images, for the viewer to load and display. Displaying cave surveys, which are zoomed out, is a problem that was hard to solve. They either become pixelated, or become very faint due to aliasing. This was solved, by downscaling the original bitmap, taking the maximum pixel value to all the prescaled pixels to make the zoomed out pixel. This obviously produces a very blocky image, but after suitable rotation and zooming out by a futher factor of four using antialiasing, the result seems pretty good.
| | |
| Anti-aliased | Max pixel value | Max pixel value and anti-aliased |
Rotating and scaling images to the tile, is a rather costly operation, whilst merging tiles together is relatively cheap. Thus each tile is made individually for each cave, giving the flexability to easily add and remove cave surveys. The making of the tiles is performed lazily, hence the slowness, that can be observed if you are the first person to look at a particular bit of cave at a particular zoom level.
#!/usr/bin/python
import cgi, cgitb, re
import Image, ImageDraw, cStringIO
import math
import os.path
import os
import transforms
import xml.dom.minidom
cgitb.enable()
def returnWorldTile(x, y, zoom):
dom = xml.dom.minidom.parse('cavesurveys.xml')
placements = dom.getElementsByTagName("placement")
tile = Image.new('RGBA', (256,256))
for placement in placements:
surveyBoundsFile = getSurveyDirectory(placement) + "bounds"
if not os.path.exists(surveyBoundsFile): #Lazy bounds file evaluation
img = Image.open(getSurveyDirectory(placement) + placement.attributes["image"].value)
makeBoundsFile(surveyBoundsFile, placement, img.size)
if inBounds(surveyBoundsFile, x, y, zoom):
tileFile = getSurveyDirectory(placement) + getTileFilename(x, y, zoom)
if not os.path.exists(tileFile): #Lazy survey file evaluation
makeSurveyTile(x, y, zoom, placement)
img = Image.open(tileFile)
tile.paste(img, None, img) # Merge all visible surveys together
#write to file object
f = cStringIO.StringIO()
tile.save(f, "PNG")
f.seek(0)
#output to browser
print "Content-type: image/png
"
print f.read()
def makeBoundsFile(surveyBoundsFile, placement, size):
sD = getSurveyDirectory(placement)
boundsFile = open(surveyBoundsFile, "w")
boundsFile.write(str(transforms.getImageLimits(placement, size)))
boundsFile.close()
def inBounds(surveyBoundsFile, x, y, zoom):
boundsFile = open(surveyBoundsFile, "r")
limits = eval(boundsFile.readline())
tileRect = transforms.getTileRect(x, y, zoom)
if ((tileRect[1] < limits[3]) and
(tileRect[3] > limits[2]) and
(tileRect[0] < limits[1]) and
(tileRect[2] > limits[0])):
return True
def makeSurveyTile(x, y, zoom, placement):
(nwx,nwy,sex,sey,swx,swy,nex,ney) = transforms.getTileCorners(placement, x, y, zoom)
len = math.sqrt((nwx-swx)**2 + (nwy- swy)**2)
dist = int(math.log(len, 2)) + 1
if dist <= 8:
img = Image.open(getSurveyDirectory(placement) + placement.attributes["image"].value)
transimg = img.transform((256,256), Image.QUAD, (nwx,nwy,swx,swy,sex,sey,nex,ney), Image.BICUBIC)
else:
if dist <=10:
img = Image.open(getSurveyDirectory(placement) + placement.attributes["image"].value)
transimg = img.transform((2**dist,2**dist), Image.QUAD, (nwx,nwy,swx,swy,sex,sey,nex,ney))
else:
scale = 2 ** (dist - 10)
filename = getSurveyDirectory(placement) + str(scale) + placement.attributes["image"].value
if os.path.exists(filename):
img = Image.open(filename)
transimg = img.transform((1024,1024), Image.QUAD, (nwx/scale,nwy/scale,swx/scale,swy/scale,sex/scale,sey/scale,nex/scale,ney/scale))
else:
img = Image.open(getSurveyDirectory(placement) + placement.attributes["image"].value)
transimg = img.transform((1024,1024), Image.QUAD, (nwx,nwy,swx,swy,sex,sey,nex,ney))
transimg = transimg.resize((256, 256), Image.ANTIALIAS)
transimg.save(getSurveyDirectory(placement) + getTileFilename(x, y, zoom), optimized=True)
def getTileFilename(x, y, zoom):
return 'x' + str(x) + 'y' + str(y) + 'z' + str(zoom) + '.png'
def getSurveyDirectory(placement):
return "/home/mjg/public_html/surveys/" + placement.attributes["image"].value + "/"
if __name__ == "__main__":
fs = cgi.FieldStorage()
if 'xy' in fs.keys() and 'zoom' in fs.keys():
data = re.match('\((.*), (.*)\)', fs['xy'].value)
x = int(data.group(1))
y = int(data.group(2))
zoom = fs['zoom'].value
returnWorldTile(x, y, zoom)
To line up a cave survey, at least one point on the survey needs to be known, along with the surveys scale and orientation.
The other significant side to the code is the coordinate transforms. Points such as entrances are put into google using WGS84 LatLong coordinates. However the tiles are implemented using a transverse mercator system base on pixels. To add to the confusion, cave surveys are typically in a seperate coordinate system. And another coordinate transform in required to get get from the coordinate sytem the cave survey is in to pixel values of the bitmap.
Some of these coordinate transforms use the GDAL library, other algorithems were seperately written or found lying around the web. At the moment transverse merctor projections can be specified and uploaded on to the site. If you require other coordinate systems, it may be possible to implement them.
#!/usr/bin/python
import osr
import re
import math
import xml.dom.minidom
import copy
import cgi
import pickle
osletters = {'A':(0,4), 'B':(1,4), 'C':(2,4), 'D':(3,4), 'E':(4,4),
'F':(0,3), 'G':(1,3), 'H':(2,3), 'J':(3,3), 'K':(4,3),
'L':(0,2), 'M':(1,2), 'N':(2,2), 'O':(3,2), 'P':(4,2),
'Q':(0,1), 'R':(1,1), 'S':(2,1), 'T':(3,1), 'U':(4,1),
'V':(0,0), 'W':(1,0), 'X':(2,0), 'Y':(3,0), 'Z':(4,0)}
def OSGridLetters2NE(GridRef):
for precision in range(1,10):
p = str(precision)
match = re.search('\s*([a-zA-Z])\s*([a-zA-Z])\s*(\d{' + p + ',' + p + '})\s*(\d{' + p + ',' + p + '})\s*\Z', GridRef)
if match:
return((osletters[match.groups()[0]][0] - 2) * 500000 + osletters[match.groups()[1]][0] * 100000 + int(match.groups()[2]) * 10 ** (5 - precision),
(osletters[match.groups()[0]][1] - 1) * 500000 + osletters[match.groups()[1]][1] * 100000 + int(match.groups()[3]) * 10 ** (5 - precision))
def getSr(name):
datumfile = open("../datums/" + name, 'r')
fs = pickle.load(datumfile)
datumfile.close()
sr = osr.SpatialReference()
if fs['datum'].value == 'other':
pmoffset = 0
if fs['primeMeridian'].value == 'Other':
pmoffset = float(fs['PMoffset'].value)
sr.SetGeogCS( "My coordinate system", "My Datum", "My Ellipsoid",
float(fs['SemiMajor'].value), float(fs['InvFlatterning'].value),
fs['primeMeridian'].value, pmoffset,
"degree" , 0.0174532925199433);
sr.SetTOWGS84(float(fs['Xoffset'].value), float(fs['Yoffset'].value), float(fs['Zoffset'].value),
float(fs['Xrot'].value), float(fs['Yrot'].value), float(fs['Zrot'].value),
float(fs['scaling'].value))
else:
sr.SetWellKnownGeogCS(fs['datum'].value);
if fs['grid'].value == 'LatLng':
pass
elif fs['grid'].value == 'UTM':
sr.SetUTM( int(fs['zone'].value), True );
elif fs['grid'].value == 'UTMsouth':
sr.SetUTM( int(fs['zone'].value), False );
elif fs['grid'].value == 'TM':
sr.SetTM(float(fs['clat'].value), float(fs['clng'].value), float(fs['scale'].value),
float(fs['falseeasting'].value), float(fs['falsenorthing'].value))
if fs['format'].value == 'osgrid':
pass
return sr
wgs = osr.SpatialReference()
wgs.SetProjCS( "WGS84" )
wgs.SetWellKnownGeogCS("WGS84")
def getTransToWGS84(name):
spatialReference = getSr(name)
return osr.CoordinateTransformation(spatialReference, wgs)
def getTransFromWGS84(name):
spatialReference = getSr(name)
return osr.CoordinateTransformation(wgs, spatialReference)
def getTileRect(x, y, zoom):
zoom = int(zoom)
x = int(x)
y = int(y)
tilesAtThisZoom = 2 ** zoom
lngWidth = 360.0 / tilesAtThisZoom # width in degrees longitude
lng = -180 + (x * lngWidth) # left edge in degrees longitude
latHeightMerc = 1. / tilesAtThisZoom # height in "normalized" mercator 0,0 top left
topLatMerc = y * latHeightMerc # top edge in "normalized" mercator 0,0 top left
bottomLatMerc = topLatMerc + latHeightMerc
# convert top and bottom lat in mercator to degrees
# note that in fact the coordinates go from about -85 to +85 not -90 to 90!
bottomLat = (180 / math.pi) * ((2 * math.atan(math.exp(math.pi * (1 - (2 * bottomLatMerc)))))
- (math.pi / 2));
topLat = (180 / math.pi) * ((2 * math.atan(math.exp(math.pi * (1 - (2 * topLatMerc))))) - (math.pi / 2))
return (lng, bottomLat, lng + lngWidth, topLat)
def getTileCorners(node, x, y, zoom):
imageSpatialReference = node.attributes["spatialreference"].value
pointnodes = node.getElementsByTagName("point")
centrePoint = readNode(pointnodes[0])
nspointnodes = node.getElementsByTagName("nspoint")
nsPoint = readNode(nspointnodes[0])
ewpointnodes = node.getElementsByTagName("ewpoint")
ewPoint = readNode(ewpointnodes[0])
(westLng, southLat, eastLng, northLat) = getTileRect(x, y, zoom)
trans = getTransFromWGS84(imageSpatialReference)
nwsr = trans.TransformPoint(westLng*math.pi/180, northLat*math.pi/180)
sesr = trans.TransformPoint(eastLng*math.pi/180, southLat*math.pi/180)
nesr = trans.TransformPoint(eastLng*math.pi/180, northLat*math.pi/180)
swsr = trans.TransformPoint(westLng*math.pi/180, southLat*math.pi/180)
nwPixels = transformToImage(nwsr, centrePoint, nsPoint, ewPoint)
sePixels = transformToImage(sesr, centrePoint, nsPoint, ewPoint)
swPixels = transformToImage(swsr, centrePoint, nsPoint, ewPoint)
nePixels = transformToImage(nesr, centrePoint, nsPoint, ewPoint)
return(nwPixels[0], nwPixels[1], sePixels[0], sePixels[1], swPixels[0], swPixels[1], nePixels[0], nePixels[1])
def getImageLimits(node, size):
imageSpatialReference = node.attributes["spatialreference"].value
pointnodes = node.getElementsByTagName("point")
centrePoint = readNode(pointnodes[0])
nspointnodes = node.getElementsByTagName("nspoint")
nsPoint = readNode(nspointnodes[0])
ewpointnodes = node.getElementsByTagName("ewpoint")
ewPoint = readNode(ewpointnodes[0])
trans = getTransToWGS84(imageSpatialReference)
bl = transformFromImage((0, 0), centrePoint, nsPoint, ewPoint)
br = transformFromImage((0, size[1]), centrePoint, nsPoint, ewPoint)
ul = transformFromImage((size[0], 0), centrePoint, nsPoint, ewPoint)
ur = transformFromImage((size[0], size[1]), centrePoint, nsPoint, ewPoint)
bl = trans.TransformPoint(bl[0], bl[1])
br = trans.TransformPoint(br[0], br[1])
ul = trans.TransformPoint(ul[0], ul[1])
ur = trans.TransformPoint(ur[0], ur[1])
return (min(bl[0], br[0], ul[0], ur[0]) * 180 / math.pi,
max(bl[0], br[0], ul[0], ur[0]) * 180 / math.pi,
min(bl[1], br[1], ul[1], ur[1]) * 180 / math.pi,
max(bl[1], br[1], ul[1], ur[1]) * 180 / math.pi)
def readNode(pointnode):
dic = {}
dic['ImX'] = float(pointnode.attributes["imx"].value)
dic['ImY'] = float(pointnode.attributes["imy"].value)
dic['GridNorth'] = float(pointnode.attributes["gridnorth"].value)
dic['GridEast'] = float(pointnode.attributes["grideast"].value)
return dic
def transformToImage(sr, centrePoint, nsPoint, ewPoint):
assert centrePoint['GridEast'] == nsPoint['GridEast']
assert centrePoint['GridNorth'] == ewPoint['GridNorth']
northSr = centrePoint['GridNorth'] - nsPoint['GridNorth']
eastSr = centrePoint['GridEast'] - ewPoint['GridEast']
northPixVec = ((centrePoint['ImX'] - nsPoint['ImX'])/northSr,
(centrePoint['ImY'] - nsPoint['ImY'])/northSr)
eastPixVec = ((centrePoint['ImX'] - ewPoint['ImX'])/eastSr,
(centrePoint['ImY'] - ewPoint['ImY'])/eastSr)
offsetSr = (sr[0]-centrePoint['GridEast'], sr[1]-centrePoint['GridNorth'])
offsetPixels = (offsetSr[0] * eastPixVec[0] + offsetSr[1] * northPixVec[0],
offsetSr[0] * eastPixVec[1] + offsetSr[1] * northPixVec[1])
return (centrePoint['ImX'] + offsetPixels[0], centrePoint['ImY'] + offsetPixels[1])
def transformFromImage(xy, centrePoint, nsPoint, ewPoint):
assert centrePoint['GridEast'] == nsPoint['GridEast']
assert centrePoint['GridNorth'] == ewPoint['GridNorth']
northSr = centrePoint['GridNorth'] - nsPoint['GridNorth']
eastSr = centrePoint['GridEast'] - ewPoint['GridEast']
eastPixVec = ((centrePoint['ImX'] - ewPoint['ImX'])/eastSr,
(centrePoint['ImY'] - ewPoint['ImY'])/eastSr)
northPixVec = ((centrePoint['ImX'] - nsPoint['ImX'])/northSr,
(centrePoint['ImY'] - nsPoint['ImY'])/northSr)
determinant = eastPixVec[0] * northPixVec[1] - eastPixVec[1] * northPixVec[0]
eSrVec = ( northPixVec[1] / determinant, -eastPixVec[1] / determinant )
nSrVec = ( -northPixVec[0] / determinant, eastPixVec[0] / determinant )
offsetPixels = ( xy[0] - centrePoint['ImX'], xy[1] - centrePoint['ImY'] )
offsetSr = (offsetPixels[0] * eSrVec[0] + offsetPixels[1] * nSrVec[0],
offsetPixels[0] * eSrVec[1] + offsetPixels[1] * nSrVec[1])
return (centrePoint['GridEast'] + offsetSr[0], centrePoint['GridNorth'] + offsetSr[1])