Message Boards Message Boards

GROUPS:

Wolfram Language OSM map extract

Posted 4 years ago
9681 Views
|
6 Replies
|
24 Total Likes
|

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.

6 Replies

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

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 4 years ago

Ha - very cool!

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

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:

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

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