13
|
16662 Views
|
6 Replies
|
25 Total Likes
View groups...
Share
GROUPS:

# Wolfram Language OSM map extract

Posted 9 years ago
 I've been trying to procedurally generate schematics of towns and cities. I turned to Mathematica to extract exemplar data to train my algorithms. I thought I'd share the data extraction process. Open street map has some really interesting data. For example, this area of Kunduz Afghanistan has high quality building polygon data, exactly what I was after. One can only conjecture why high density LIDAR data of this area exists, or how it ended up so accessible. OSM has an export button, provided you select a small enough region, it saves you an '.osm' file. Two example files from Kunduz are attached. .osm primarily stores two structures: nodes (single geographic points) and ways (lists of nodes that form polygons, tagged to reflect content type). We can easily drag these out in Mathematica to sift out our target data. First importing a data set, cutting an import string into individual sections. Then selecting the node type data, and saving the node coordinates as functions of their IDs, and saving those IDs in 'nodes'. strings = StringDrop[StringTrim[#], 1] & /@ StringSplit[Import["C:\\Users\\Me\\Downloads\\map.osm"], ">"]; nodestrings = Cases[If[StringSplit[#][[1]] == "node", StringDrop[#, 5]] & /@ strings, Except[Null]]; (node[#[[1]]] = {#[[3]], #[[2]]}; AppendTo[nodes, #[[1]]];) & /@ (ToExpression[StringDrop[StringTrim[StringDrop[#[[1]], #[[2]]]], -#[[3]]] & /@ Transpose[{StringSplit[#][[{1, -2, -1}]], {4, 5, 5}, {1, 1, 2}}] & /@ nodestrings]);)  Next I select the non node data, extracting ways, their node references, and their tags; stored in 'waylists' and 'taglists'. nonnode = Cases[If[StringSplit[#][[1]] != "node", #] & /@ (StringDrop[StringTrim[#], 1] & /@ StringSplit[Import["C:\\Users\\Me\\Downloads\\map.osm"], ">"]), Except[Null]]; openers = Parallelize[StringSplit[#][[1]] & /@ # & /@ SplitBy[nonnode, StringSplit[#][[1]] &]]; nonnode = SplitBy[nonnode, StringSplit[#][[1]] &]; AppendTo[waylists, Parallelize[nonnode[[#[[1]] ;; #[[2]]]] & /@ Transpose[{First /@ Position[openers, "way"], First /@ Position[openers, "/way"]}]]]; AppendTo[taglists, Cases[If[StringTake[#, 1] == "k", StringDrop[StringDrop[#, 3], -1]] & /@ #, Except[Null]][[1]] & /@ StringSplit[#[[3]]] & /@ waylists[[-1]]];  Now the data can be easily polled. For example, here all ways with the building tag are dumped into lists of locations, allowing easy visualisation of building polygons. buildpolys = Table[(node[#] & /@ ToExpression[StringDrop[StringDrop[#, 8], -2] & /@ #[[2]]]) & /@ waylists[[i]][[First /@ Position[If[MemberQ[#, "building"], 1] & /@ taglists[[i]], 1]]], {i, 1, 1}]; Graphics[{Blue, Line[#] & /@ buildpolys}]  All these functions were designed to handle multiple files, allowing larger imports. See the attached notebook for an example. Attachments:
6 Replies
Sort By:
Posted 9 years ago
 Very interesting! Do you think it is possible to overlay the OSMs somehow on GeoGraphics? Are they in the same geo-coordinate system?
Posted 9 years ago
 Using GeoGraphics ?! Of course - great idea! Things then become even easier: The geographic coordinates (latitude/longitude) can be used directly as extracted from the osm-file, no further transformation is necessary. Some minimal code: geoBuildingLines = << "geoBuildingLines.txt"; GeoGraphics[GeoPath[geoBuildingLines]] Nice things are possible (here just an appetizer):Henrik Attachments:
Posted 9 years ago
 OSM provides a handful of APIs. You can easily summon the relevant file in code, allowing work entirely in Wolfram Geo functions. (* Establish area bounds *) GeoLoc = FindGeoLocation[]; ValLoc = GeoLoc[[1]]; MetersRadius = MR = 1000; Corners = ToString /@ Join[ GeoDestination[GeoDestination[GeoLoc, {MR, 180}], {MR, 270}][[1]][[{2, 1}]], GeoDestination[GeoDestination[GeoLoc, {MR, 0}], {MR, 90}][[1]][[{2, 1}]]]; (* request relevent map *) ImportRegion = URLFetch[StringJoin[ Flatten[{"http://api.openstreetmap.org/api/0.6/map?bbox=", Transpose[{Corners, {",", ",", ",", ""}}]}]]]; (* previously seen unpackaging *) nodes = {}; Importnodes[ImportRegion]; waylists = {}; taglists = {}; Importways[ImportRegion]; (* extract desired data *) buildpolys = Quiet[Table[(node[#] & /@ ToExpression[StringDrop[StringDrop[#, 8], -2] & /@ #[[2]]]) & /@ waylists[[i]][[First /@ Position[If[MemberQ[#, "building"], 1] & /@ taglists[[i]],1]]], {i, 1, 2}]]; (* plot onto map *) GeoGraphics[{GeoPath[Reverse /@ # & /@ buildpolys[[1]], "Rhumb"], GeoMarker[]}] 
Posted 9 years ago
 Ha - very cool!
Posted 9 years ago
 Dear David,thank you very much for sharing this really nice idea! I already had a lot of fun playing around with this new possibility! I would like to make some minor remarks: I found the coding much easier when the osm-file is imported as "XML" (as it is) instead of "TXT"; a useful addition to your code above could be a transformation from geographic coordinates (latitude, longitude) to planimetric coordinates (e.g. UTM grid). For this I am using the following little functions: . (* get UTMZone from geographic longitude *) getUTM[lon_] := Module[{x, index}, x = If[lon > 180, lon - 360, lon]; index = Round[(x + 3.)/6.] + 30; StringJoin["UTMZone", ToString[index]]] (* calculate: grad to UTM *) grad2utm[latlong_List, utmZone_String] := First@GeoGridPosition[GeoPosition[latlong, "WGS84"], GeoProjectionData[utmZone]] (* calculate: UTM to grad *) utm2grad[xy_List, utmZone_String] := QuantityMagnitude@ LatitudeLongitude[GeoGridPosition[xy, GeoProjectionData[utmZone]]] Regards Henrik
Posted 9 years ago
 Thanks for pointing out the XML basis of OSM files. I probably burnt a lot of time overall by not looking the format up properly, just hacking through it.