Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
592 lines (417 sloc) 14.1 KB
#! /usr/local/bin/python
import voronoi
import sys
import getopt
from pylab import *
#import matplotlib
from shapely.geometry import Polygon
from shapely.ops import polygonize
import hashlib
import json
import time
import globalmaptiles
WorldRanges={}
WorldRanges["AUSTIN"]=[30.8, -98.5, 29.535, -97.031]
WorldRanges["TX"]=[36.5, -106, 25, -93]
WorldRanges["US"]=[55, -130, 23, -60]
WorldRanges["GUS"]=[60, -140, 22, -50]
WorldRanges["KR"]=[45, 120, 32, 135]
WorldRanges["W"]=[90, -180, -90, 180]
World="W"
WorldRange=WorldRanges[World]
PlotIt=False
def checkInRange(WorldRange,x,y):
if x>=WorldRange[1] and x<=WorldRange[3]:
if y<=WorldRange[0] and y>=WorldRange[2]:
return True
return False
def getExtremes(line, m_range):
a,b,c=line
if a==0:
y0=(c-(a*m_range["min_x"]))/b
y1=(c-(a*m_range["max_x"]))/b
return [(m_range["min_x"],y0), (m_range["max_x"], y1)]
if b==0:
x0=(c-(b*m_range["min_y"]))/a
x1=(c-(b*m_range["max_y"]))/a
return [(x0, m_range["min_y"]), (x1, m_range["max_y"])]
x0=(c-(b*m_range["min_y"]))/a
x1=(c-(b*m_range["max_y"]))/a
y0=(c-(a*m_range["min_x"]))/b
y1=(c-(a*m_range["max_x"]))/b
return [(x0, m_range["min_y"]), (x1, m_range["max_y"]), (m_range["min_x"],y0), (m_range["max_x"], y1)]
def getExtreme(line, v_known, LR):
global WorldRange
global_extremes=getExtremes(line, {"min_y":WorldRange[2], "max_y":WorldRange[0], "min_x":WorldRange[1], "max_x":WorldRange[3]})
v_unknown=None
CurrentExtreme=None
for extreme in global_extremes:
if checkInRange(WorldRange,extreme[0],extreme[1])==False:
#v_unknown=CurrentExtreme
continue
if v_known[0]<extreme[0] and LR==1:
v_unknown=extreme
if v_known[0]>extreme[0] and LR==0:
v_unknown=extreme
if CurrentExtreme:
if CurrentExtreme[0]==extreme[0]:
if CurrentExtreme:
if CurrentExtreme[1]<extreme[1] and LR==0:
v_unknown=extreme
else:
v_unknown=CurrentExtreme
CurrentExtreme=extreme
#print LR, v_known, extreme, v_unknown
return v_unknown
def update_maxmin(m_range, x, y):
if (m_range["max_x"]<x):
m_range["max_x"]=x
if (m_range["max_y"]<y):
m_range["max_y"]=y
if (m_range["min_x"]>x):
m_range["min_x"]=x
if (m_range["min_y"]>y):
m_range["min_y"]=y
return m_range
def linkExtremes(point1, point2, m_range):
#print "Link Extremes"
#print point1, point2
global WorldRange
if point1[1]==point2[1] or point1[0]==point2[0]:
#print "already done case"
return [(point1, point2)]
if point1[1]==-point2[1] or point1[0]==-point2[0]:
#print "worst case"
if point1[1]==-point2[1]:
if abs(point1[0]-m_range["min_x"])+abs(point2[0]-m_range["min_x"])<abs(point1[0]-m_range["max_x"])+abs(point2[0]-m_range["max_x"]):
output=[(point1, (WorldRange[1], point1[1])), ((WorldRange[1], point1[1]),(WorldRange[1], point2[1])), ((WorldRange[1], point2[1]), point2)]
else:
output=[(point1, (WorldRange[3], point1[1])), ((WorldRange[3], point1[1]),(WorldRange[3], point2[1])), ((WorldRange[3], point2[1]), point2)]
else:
if abs(point1[1]-m_range["min_y"])+abs(point2[1]-m_range["min_y"])<abs(point1[1]-m_range["max_y"])+abs(point2[1]-m_range["max_y"]):
output=[(point1, (point1[0], WorldRange[2])), ((point1[0], WorldRange[2]),(point2[0], WorldRange[2])), ((point2[0], WorldRange[2]), point2)]
else:
output=[(point1, (point1[0], WorldRange[0])), ((point1[0], WorldRange[0]),(point2[0], WorldRange[0])), ((point2[0], WorldRange[0]), point2)]
return output
#print "one corner case"
if point1[1]==WorldRange[0] or point1[1]==WorldRange[2]:
#print "add x:", point2[0]
#print "add y:", point1[1]
output=[(point1, (point2[0], point1[1])), ((point2[0], point1[1]), point2)]
else:
#print "add x:", point1[0]
#print "add y:", point2[1]
output=[(point1, (point1[0], point2[1])), ((point1[0], point2[1]), point2)]
return output
def VoronoiLineEdges(PointsMap):
Sitepts = []
pts = {}
#print CurrentDate, PointsMap[PointsMap.keys()[0]]
for grid, stn in PointsMap.items():
x=float(stn[0])
y=float(stn[1])
station=grid
#station.extend( stn[3:])
#print x,y,station
pts[ (x,y) ]=station
stncounts=len(pts.keys())
#print stncounts, "points"
site_points=[]
for pt in pts.keys():
Sitepts.append(voronoi.Site(pt[0],pt[1]))
site_points.append( (pt[0],pt[1]) )
#print "Calculating Voronoi Lattice",
siteList = voronoi.SiteList(Sitepts)
context = voronoi.Context()
voronoi.Edge.EDGE_NUM=0
voronoi.voronoi(siteList,context)
vertices=context.vertices
lines=context.lines
edges=context.edges
triangles=context.triangles
has_edge=context.has_edge
return vertices, lines, edges, has_edge
def VoronoiGeoJson_MultiPolygons(PointsMap, BoundingBox="W",PlotMap=False):
vl=VoronoiPolygons(PointsMap, BoundingBox="W", PlotMap=PlotMap)
cluster_id="VoronoiLattice"
geojson={}
geojson["properties"]={}
geojson["properties"]["_domain_id"]=cluster_id
geojson["type"]="Feature"
m = hashlib.md5()
hash_string=cluster_id
m.update(hash_string)
geojson["id"]=m.hexdigest()
geojson["geometry"]={}
geojson["geometry"]["type"]="MultiPolygon"
geojson["geometry"]["coordinates"]=[]
for vl_num, data in vl.items():
geojson["geometry"]["coordinates"].append([list(data["obj_polygon"].exterior.coords)])
return json.dumps(geojson)
def VoronoiGeoJson_Polygons(PointsMap, BoundingBox="W",PlotMap=False):
vl=VoronoiPolygons(PointsMap, BoundingBox="W", PlotMap=PlotMap)
output=""
for vl_num, data in vl.items():
#print data
cluster_id=data["info"]
geojson={}
geojson["properties"]={}
geojson["properties"]["_domain_id"]=cluster_id
geojson["type"]="Feature"
m = hashlib.md5()
hash_string=cluster_id
m.update(hash_string)
geojson["id"]=m.hexdigest()
geojson["geometry"]={}
geojson["geometry"]["type"]="Polygon"
coords=[]
geojson["geometry"]["coordinates"]=[list(data["obj_polygon"].exterior.coords)]
output+=json.dumps(geojson)#, sort_keys=True, indent=3)
output+="\n"
return output
def VoronoiPolygons(PointsMap, BoundingBox="W", PlotMap=False):
global PlotIt
global WorldRange
global WorldRanges
if type(BoundingBox)==type([]):
if len(BoundingBox)==4:
WorldRange=BoundingBox
else:
return "Error in Bounding Box"
else:
WorldRange=WorldRanges[BoundingBox]
PlotIt=PlotMap
currenttime=time.time()
Sitepts = []
pts = {}
#print CurrentDate, PointsMap[PointsMap.keys()[0]]
for grid, stn in PointsMap.items():
x=float(stn[0])
y=float(stn[1])
station=grid
#station.extend( stn[3:])
#print x,y,station
pts[ (x,y) ]=station
stncounts=len(pts.keys())
#print stncounts, "points"
site_points=[]
for pt in pts.keys():
Sitepts.append(voronoi.Site(pt[0],pt[1]))
site_points.append( (pt[0],pt[1]) )
#print "Calculating Voronoi Lattice",
siteList = voronoi.SiteList(Sitepts)
context = voronoi.Context()
voronoi.Edge.EDGE_NUM=0
voronoi.voronoi(siteList,context)
vertices=context.vertices
lines=context.lines
edges=context.edges
#print edges
#For Faster Access
edge_dic={}
for edge in edges:
edge_dic[edge[0]]=edge[1:]
triangles=context.triangles
has_edge=context.has_edge
voronoi_lattice={}
m_range={}
m_range["max_x"]=-9999999999
m_range["max_y"]=-9999999999
m_range["min_x"]=9999999999
m_range["min_y"]=9999999999
#Get the range!!
for pnt in site_points:
m_range=update_maxmin(m_range, pnt[0], pnt[1])
#print "Getting the Polygons"
prev_percent=0
for station, ls in has_edge.items():
voronoi_lattice[station]={}
voronoi_lattice[station]["coordinate"]=site_points[station]
voronoi_lattice[station]["info"]=pts[ site_points[station] ]
polygon=[]
prev_extreme=[]
Verbose=True
if Verbose:
current_percent=int(station/float(stncounts)*100)
if current_percent!=prev_percent:
#print station,"/", stncounts, current_percent, "% Done"
timeelapse=time.time()-currenttime
#print station, timeelapse
currenttime=time.time()
prev_percent=current_percent
#For every lines that the station owns
for l in ls:
e=edge_dic[l]
v1=vertices[e[0]]
v2=vertices[e[1]]
if e[0] < 0 and checkInRange(WorldRange,v2[0],v2[1])==False: continue
if e[1] < 0 and checkInRange(WorldRange,v1[0],v1[1])==False: continue
if e[0] > -1 and e[1] > -1 and checkInRange(WorldRange,v1[0],v1[1])==False and checkInRange(WorldRange,v2[0],v2[1])==False:
continue
if e[0] < 0 or checkInRange(WorldRange,v1[0],v1[1])==False:
v1=getExtreme(lines[l],v2, LR=0)
if len(prev_extreme)==0:
prev_extreme=v1
else:
extreme_points=linkExtremes(prev_extreme, v1, m_range)
for extreme_pair in extreme_points:
polygon.append(extreme_pair)
if e[1] < 0 or checkInRange(WorldRange,v2[0],v2[1])==False :
v2=getExtreme(lines[l],v1, LR=1)
if len(prev_extreme)==0:
prev_extreme=v2
else:
extreme_points=linkExtremes(prev_extreme, v2, m_range)
for extreme_pair in extreme_points:
polygon.append(extreme_pair)
if v1!=v2:
polygon.append( (v1,v2) )
if len(polygon)==0:
sys.stderr.write ("\nThis station does not have meaningful polygon:")
sys.stderr.write (str(pts[ site_points[station] ])+" at "+str(site_points[station])+"\n")
for l in ls:
e=edge_dic[l]
v1=vertices[e[0]]
v2=vertices[e[1]]
sys.stderr.write(str(e[0])+","+str(e[1])+"\n")
sys.stderr.write(str(v1)+","+str(v2)+"\n")
voronoi_lattice.pop(station)
continue
#print polygon
try:
result = list(polygonize( polygon ))
except ValueError:
#print "Wrong:", polygon
#Delete invisible short lines
point_to_replace=[]
for line in polygon:
d=math.hypot( line[0][0]-line[1][0], line[0][1]-line[1][1])
if d<1:
polygon.remove(line)
point_to_replace=line
i=0
New_Lines=[]
for line in polygon:
New_Lines.append(list(line))
j=0
for point in line:
if point==point_to_replace[0]:
#print line, point_to_replace
New_Lines[i][j]=point_to_replace[1]
j+=1
i+=1
polygon=tuple(New_Lines)
#I could not figure out why sometimes it fails to draw a polygon.
try:
result = list(polygonize( polygon ))
except:
voronoi_lattice.pop(station)
continue
finalpoly=result[0]
#print list(finalpoly.exterior.coords)
#voronoi_lattice[station]["polygon"]=list(finalpoly.exterior.coords)
voronoi_lattice[station]["obj_polygon"]=finalpoly
#print polygon
if (PlotIt):
for station, data in voronoi_lattice.items():
x=[]
y=[]
#print data
try:
polygon_data=data["obj_polygon"]
#print polygon_data
for point in list(polygon_data.exterior.coords):
x.append(point[0])
y.append(point[1])
except:
print "Error", data["name"]
#if station == 8 :
#plot(x,y)
fill(x,y, alpha=0.6)
plot(data["coordinate"][0],data["coordinate"][1])
#text(data["coordinate"][0],data["coordinate"][1],data["info"])
show()
return voronoi_lattice
ClosestZoomLevel=7
def plot_voronoi(voronoi_lattice):
for station, data in voronoi_lattice.items():
x=[]
y=[]
polygon_data=data["obj_polygon"]
#print polygon_data
for point in list(polygon_data.exterior.coords):
x.append(point[0])
y.append(point[1])
#if station == 8 :
#plot(x,y)
fill(x,y, alpha=0.3)
#scatter(data["coordinate"][0], data["coordinate"][1])
# show()
def quadGrid(grid_range):
#print grid_range
bix=(grid_range[0]+grid_range[2])/2
biy=(grid_range[1]+grid_range[3])/2
grid0=[grid_range[0], grid_range[1], bix, biy]
grid1=[grid_range[0], biy, bix, grid_range[3]]
grid2=[bix, grid_range[1], grid_range[2], biy]
grid3=[bix, biy, grid_range[2], grid_range[3]]
return (grid0, grid1, grid2, grid3)
def polygonize_grid(grid_range):
x0=grid_range[1]
y0=grid_range[0]
x1=grid_range[3]
y1=grid_range[2]
polygon=( (x0, y0), (x0, y1), (x1, y1), (x1, y0), (x0, y0))
return polygon
def get_quadkeystr(quadkey):
for i in range(ClosestZoomLevel-len(quadkey)):
quadkey+="X"
return quadkey
#Zoom in!!
def GridMap(quadkey="", ploygons_to_lookat=None):
mercator = globalmaptiles.GlobalMercator()
global PlotIt
for quadkeyadd in range(4):
curQuadKey=quadkey+str(quadkeyadd)
tx, ty, zl=mercator.QuadTree2TMS(curQuadKey)
grid_latlon=mercator.TileLatLonBounds(tx, ty, zl)
#print grid_latlon
polygon=polygonize_grid(grid_latlon)
#print polygon
obj_polygon=Polygon(polygon)
inter_cnt=0
Grid_Weather_Stations=[]
ploygons_to_lookat_new=[]
for pobj in ploygons_to_lookat:
#print pobj, obj_polygon
if pobj[0].intersection(obj_polygon):
inter_cnt+=1
ploygons_to_lookat_new.append(pobj)
#print inter_cnt
if inter_cnt> 2 and zl<ClosestZoomLevel:
GridMap(curQuadKey, ploygons_to_lookat_new)
continue
quadkeystr=get_quadkeystr(curQuadKey)
for pobj in ploygons_to_lookat_new:
if PlotIt:
x=[]
y=[]
for point in list(obj_polygon.exterior.coords):
x.append(point[0])
y.append(point[1])
plot(x,y, alpha=0.3)
print quadkeystr+"\t"+pobj[1]+"\t"+str(pobj[2][0])+"\t"+str(pobj[2][1])
def GridVoronoi(voronoi_lattice, zl, PlotMap=True):
#print "Mapping the Grids"
global PlotIt
PlotIt=PlotMap
global ClosestZoomLevel
ClosestZoomLevel=zl
if PlotMap:
plot_voronoi(voronoi_lattice)
ploygons_to_lookat_new=[]
for station, data in voronoi_lattice.items():
ploygons_to_lookat_new.append((data["obj_polygon"], data["info"], data["coordinate"] ))
GridMap("", ploygons_to_lookat_new)
if PlotMap:
show()