Message Boards Message Boards

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

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}]

Kunduz example

All these functions were designed to handle multiple files, allowing larger imports. See the attached notebook for an example.

Attachments:
POSTED BY: David Gathercole
6 Replies

Very interesting! Do you think it is possible to overlay the OSMs somehow on GeoGraphics? Are they in the same geo-coordinate system?

POSTED BY: Sam Carrettie

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):

enter image description here

Henrik

Attachments:
POSTED BY: Henrik Schachner

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[]}]

example output

POSTED BY: David Gathercole
Posted 9 years ago

Ha - very cool!

POSTED BY: Stephen Coast

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 BY: Henrik Schachner

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.

POSTED BY: David Gathercole
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract