Le 16/12/2009 12:24, Bruno Cortial a écrit :
Bonjour,
Je ne trouve pas de statistiques kilométriques sur le réseau routier
d'OSM. Du genre nb de kilomètres de residential pour la commune de
Pornic, ou de primary/secondary pour la Loire-Atlantique ? On parle de
trou dans la couverture du réseau routier, ce genre d'outil pourrait
aiguillonner les mappeurs. Emilie a déjà évoqué des outils d'avancement
OSM, cela pourrait en être une ébauche.
Il nous faudrait également une base de comparaison. Les conseils
généraux/régionaux fournissent-ils ce genre d'info ?
A titre de petite fourmi, si je devais me lancer dans ce genre de
requête, l'import PostGIS s'impose, mais je n'y connais rien: quelqu'un
a-t-il des exemples de requêtes SQL sur une base OSM importée ?
Voila un exemple n'utilisant pas de base postgis mais basé sur un
parseur sax d'un fichier osm. Il faut ensuite utiliser osmosis pour
découper un fichier osm en région/département.
Pour les requêtes postgis, je ne sait pas dans quelles projection ils
faut se placer pour avoir une distance sphérique, je laisse donc les
spécialistes te répondre.
--
Etienne
# Copyright Arno arenev...@fdn.fr
import sys, os, math
from xml.sax import make_parser, handler
import xml
# function from http://www.johndcook.com/python_longitude_latitude.html
def distance_on_unit_sphere(lat1, long1, lat2, long2):
#print "distance_on_unit_sphere %.25f %.25f %.25f %.25f" % (lat1, long1,
lat2, long2)
if (lat1 == lat2) and (long1 == long2):
return 0
# Convert latitude and longitude to
# spherical coordinates in radians.
degrees_to_radians = math.pi/180.0
# phi = 90 - latitude
phi1 = (90.0 - lat1)*degrees_to_radians
phi2 = (90.0 - lat2)*degrees_to_radians
# theta = longitude
theta1 = long1*degrees_to_radians
theta2 = long2*degrees_to_radians
# Compute spherical distance from spherical coordinates.
# For two locations in spherical coordinates
# (1, theta, phi) and (1, theta, phi)
# cosine( arc length ) =
# sin phi sin phi' cos(theta-theta') + cos phi cos phi'
# distance = rho * arc length
cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) +
math.cos(phi1)*math.cos(phi2))
arc = math.acos( cos )
# Remember to multiply arc by the radius of the earth
# in your favorite set of units to get length.
return arc
class parseOsm(handler.ContentHandler):
def __init__(self):
self.waynodes = []
self.nodes = {}
self.inaway = False
self.highwaytype = ""
self.results = {}
def loadOsm(self, filename):
if(not os.path.exists(filename)):
return
print "loading ", filename
print
try:
parser = make_parser()
parser.setContentHandler(self)
parser.parse(filename)
except xml.sax._exceptions.SAXParseException:
print "Error loading %s" % filename
def printResults(self):
items = sorted(self.results.items(), self._sortResults)
totdist = 0
for res in items:
print "%s: %.3f km" % (res[0], res[1] / 1000)
totdist += res[1] / 1000
print "--------------------------"
print "total distance: %.3f kms" % (totdist)
def _sortResults (self, tuple1, tuple2):
if (tuple1[1] < tuple2[1]):
return -1
elif (tuple1[1] > tuple2[1]):
return 1
if (tuple1[1] == tuple2[1]):
return 0
def startElement(self, name, attrs):
if name == 'way':
self.startWay()
elif name == 'node':
id = int(attrs.get('id'))
lat = float(attrs.get('lat'))
lon = float(attrs.get('lon'))
self.nodes[id] = (lat,lon)
elif name == 'nd':
if self.inaway:
ref = int(attrs.get('ref'))
self.waynodes.append(self.nodes[ref])
elif name == 'tag':
if self.inaway:
key = attrs.get('k')
value = attrs.get('v')
if key == "highway":
self.highwaytype = value
def startWay(self):
self.inaway = True
self.highwaytype = ""
self.waynodes = []
def endWay(self):
self.inaway = False
self.highwaytype = ""
self.waynodes = []
def endElement(self, name):
if name == 'way':
if len(self.highwaytype) == 0:
self.endWay()
return
if len(self.waynodes) < 2:
self.endWay()
return
dist = 0
curLat = self.waynodes[0][0]
curLon = self.waynodes[0][1]
for node in self.waynodes[1:]:
lat = node[0]
lon = node[1]
dist += distance_on_unit_sphere(curLat, curLon, lat, lon) * 6373 *
1000
curLat = lat
curLon = lon
if self.results.get(self.highwaytype):
self.results[self.highwaytype] += dist
else:
self.results[self.highwaytype] = dist
self.endWay()
if (__name__ == "__main__"):
if len(sys.argv) != 2:
print "please specify an osm file name"
sys.exit(2)
filename = sys.argv[1]
parser = parseOsm()
parser.loadOsm(filename)
parser.printResults()
_______________________________________________
Talk-fr mailing list
Talk-fr@openstreetmap.org
http://lists.openstreetmap.org/listinfo/talk-fr
