23
$\begingroup$

Mathematica is great at importing ESRI shapefiles (.shp).

Import["http://exampledata.wolfram.com/usamap.zip", "Graphics"]

enter image description here

But it seems like it's a glaring omission that it can't export them. Does MMA have the ability?

$\endgroup$
27
$\begingroup$

Yes it can. Now it's not necessarily easy, but all the groundwork is in place to build a function that can export ESRI shapefiles.

(A lot of code follows....) UPDATE: Available on GitHub as a Wolfram Language Package. See here.

First, a shapefile consists of at least three essential files: 1) the .shp main file, 2) the .shx index file, and 3) the .dbf--a hateful file format that stores the tabular data related to the feature. Sadly, the .dbf file is another one that Mathematica can import but not export. The .shp and .shx files change endianness frequently, but are otherwise pretty straight forward.

I found it easier to write functions that first import the data so I could study how it worked. I then reversed that to write the export functions. I can include those functions if there is a desire. Only the Export functions are shown here.

Now to the functions:

Here's a function to write the .shp and .shx files simultaneously:

writeshp[geometry_, filepath_] := 
 Module[{str = OpenWrite[filepath, BinaryFormat -> True], 
   shx = OpenWrite[StringReplace[filepath, ".shp" -> ".shx"], 
     BinaryFormat -> True], shapetype, bounds, recordnumber = 0},
  BinaryWrite[str, {9994, 0, 0, 0, 0, 0, filelength[geometry]}, 
   "Integer32", ByteOrdering -> 1];
  shapetype = 
   Pick[{1, 3, 5}, {Point, Line, Polygon}, 
     Commonest[geometry[[All, 0]]][[1]]][[1]];
  BinaryWrite[str, {1000, shapetype}, "Integer32", ByteOrdering -> -1];
  bounds = 
   If[shapetype == 1, MinMax /@ geometry[[All, 1]], 
    MinMax /@ Transpose[Join @@ (geometry[[All, 1]])]];
  BinaryWrite[
   str, {bounds[[1, 1]], bounds[[2, 1]], bounds[[1, 2]], 
    bounds[[2, 2]], 0., 0., 0., 0.}, "Real64", ByteOrdering -> -1];
  BinaryWrite[shx, {9994, 0, 0, 0, 0, 0, (100 + 8*Length@geometry)/2},
    "Integer32", ByteOrdering -> 1];
  BinaryWrite[shx, {1000, shapetype}, "Integer32", ByteOrdering -> -1];
  BinaryWrite[
   shx, {bounds[[1, 1]], bounds[[2, 1]], bounds[[1, 2]], 
    bounds[[2, 2]], 0., 0., 0., 0.}, "Real64", ByteOrdering -> -1];
  Which[
   shapetype == 1,
   Do[writepoint[str, shx, record, recordnumber++], {record, 
     geometry}],
   shapetype == 3,
   Do[writepolyline[str, shx, record, recordnumber++], {record, 
     geometry}],
   shapetype == 5,
   Do[writepolygon[str, shx, record, recordnumber++], {record, 
     geometry}]
   ];
  Close[str];
  Close[shx];
  ]

It has four helper functions:

writepolyline[stream_, shxstream_, linerecord_, recordnumber_] := 
 Module[{numpart = Depth@linerecord[[1]] - 2, numpoints, bounds},
  If[numpart > 1, numpoints = Total[Length /@ linerecord[[1]]]; 
   bounds = MinMax /@ Transpose[Join @@ (linerecord[[1]])], 
   numpoints = Length@linerecord[[1]]; 
   bounds = MinMax /@ Transpose[linerecord[[1]]]];
  BinaryWrite[
   shxstream, {StreamPosition[stream]/2, 
    22 + 2*numpart + 8*numpoints}, "Integer32", ByteOrdering -> 1];
  BinaryWrite[stream, {recordnumber, 22 + 2*numpart + 8*numpoints}, 
   "Integer32", ByteOrdering -> 1];
  BinaryWrite[
   stream, {3, bounds[[1, 1]], bounds[[2, 1]], bounds[[1, 2]], 
    bounds[[2, 2]], numpart, numpoints, 
    Sequence @@ Range[0, numpart - 1], 
    Sequence @@ (Flatten@linerecord[[1]])}, {"Integer32", "Real64", 
    "Real64", "Real64", "Real64", "Integer32", "Integer32", 
    Sequence @@ ConstantArray["Integer32", numpart], 
    Sequence @@ ConstantArray["Real64", numpoints*2]}, 
   ByteOrdering -> -1]
  ]

writepoint[stream_, shxstream_, pointrecord_, recordnumber_] := 
 Module[{numpart = Depth@pointrecord[[1]] - 2, numpoints, bounds},
  BinaryWrite[shxstream, {StreamPosition[stream]/2, 10}, "Integer32", 
   ByteOrdering -> 1];
  BinaryWrite[stream, {recordnumber, 10}, "Integer32", 
   ByteOrdering -> 1];
  BinaryWrite[
   stream, {1, Sequence @@ (Flatten@pointrecord[[1]])}, {"Integer32", 
    "Real64", "Real64"}, ByteOrdering -> -1]
  ]

writepolygon[stream_, shxstream_, polyrecord_, recordnumber_] := 
 Module[{numpart = Depth@polyrecord[[1]] - 2, numpoints, bounds},
  If[numpart > 1, numpoints = Total[Length /@ polyrecord[[1]]]; 
   bounds = MinMax /@ Transpose[Join @@ (polyrecord[[1]])], 
   numpoints = Length@polyrecord[[1]]; 
   bounds = MinMax /@ Transpose[polyrecord[[1]]]];
  BinaryWrite[
   shxstream, {StreamPosition[stream]/2, 
    22 + 2*numpart + 8*numpoints}, "Integer32", ByteOrdering -> 1];
  BinaryWrite[stream, {recordnumber, 22 + 2*numpart + 8*numpoints}, 
   "Integer32", ByteOrdering -> 1];
  BinaryWrite[
   stream, {3, bounds[[1, 1]], bounds[[2, 1]], bounds[[1, 2]], 
    bounds[[2, 2]], numpart, numpoints, 
    Sequence @@ Range[0, numpart - 1], 
    Sequence @@ (Flatten@polyrecord[[1]])}, {"Integer32", "Real64", 
    "Real64", "Real64", "Real64", "Integer32", "Integer32", 
    Sequence @@ ConstantArray["Integer32", numpart], 
    Sequence @@ ConstantArray["Real64", numpoints*2]}, 
   ByteOrdering -> -1]
  ]

filelength[geometry_] := 
 Module[{shapetype = 
    Pick[{1, 3, 5}, {Point, Line, Polygon}, 
      Commonest[geometry[[All, 0]]][[1]]][[1]], 
   records = Length@geometry[[All, 1]], 
   points = Length@Flatten@geometry[[All, 1]], 
   parts = Total[Depth /@ geometry[[All, 1]] - 2]},
  If[shapetype == 1, 28*records + 100, 
    52*records + 4*parts + 8*points + 100]/2
  ]

As mentioned before, the .dbf file also must be written. Here's a function to write the .dbf for 3 data types (strings, integers, and floats but it can easily be extended to others if needed):

writedbf[assoc_, filepath_] := 
 Module[{str = OpenWrite[filepath, BinaryFormat -> True], 
   vals = Values@assoc, fieldnames = Keys[assoc], fieldtypes, offsets,
    header, subrecords, recods, recordstart, records},
  fieldnames = 
   If[StringLength[#] > 10, StringTake[#, 10], #] & /@ fieldnames;
  fieldtypes = Commonest[Head /@ #][[1]] & /@ vals;
  Do[If[fieldtypes[[i]] == Real, 
    vals[[i]] = realformat[vals[[i]]]], {i, Length@vals}];
  offsets = 
   Table[Switch[fieldtypes[[i]], 
     String, {Max[50, Max[Length /@ vals[[i]]]], 0}, Integer, {16, 0},
      Real, {19, 11}], {i, Length@vals}];
  fieldtypes = 
   Pick[{"C", "N", "F"}, {String, Integer, Real}, #][[1]] & /@ 
    fieldtypes;
  subrecords = 
   Flatten@Table[
     Flatten[{PadRight[ToCharacterCode[fieldnames[[i]]], 11, 0], 
       ToCharacterCode[fieldtypes[[i]]], 0, 0, 0, 0, 
       Sequence @@ offsets[[i]], 0, ConstantArray[0, 13]}], {i, 
      Length@fieldnames}];
  recordstart = 32*(Length@vals + 1) + 1;
  records = 
   Transpose@
    Table[If[fieldtypes[[i]] == "N", PadLeft, PadRight][
      ToCharacterCode[ToString[vals[[i, j]]]], offsets[[i, 1]], 
      32], {i, Length@vals}, {j, Length@vals[[i]]}];
  header = {3, DateList[][[1]] - 1900, DateList[][[2]], 
    DateList[][[3]], Length@vals[[1]], recordstart, 
    Length@Flatten@records[[1]] + 1, Sequence @@ ConstantArray[0, 17],
     87, 0, 0};
  BinaryWrite[str, 
   header, {Sequence @@ ConstantArray["Byte", 4], "Integer32", 
    "Integer16", "Integer16", Sequence @@ ConstantArray["Byte", 20]}];
  BinaryWrite[str, subrecords, "Byte"];
  BinaryWrite[str, {13}, "Byte"];
  BinaryWrite[str, Flatten[Prepend[#, 32] & /@ (Join @@@ records)], 
   "Byte"];
  BinaryWrite[str, {26}, "Byte"];
  Close[str];
  ]

And one helper function to format Real numbers:

realformat[num_] := 
 ToString@ScientificForm@PaddedForm[num, {12, 11}, 
  NumberFormat -> (Row[{#1, "e", If[ToExpression[(#3 /. "" -> "0")] < 0, "-", "+"], StringPadLeft[StringReplace[ToString@#3, "-" -> ""], 3,  "0"]}] &)]
SetAttributes[realformat, Listable]

And finally a function to wrap it all together:

exportshapefile[filepath_, geometry_, assoc_] := Module[{},
  If[! StringMatchQ[filepath, __ ~~ "SHP", IgnoreCase -> True], 
   Abort[]];
  If[Length@First@Values@assoc != Length@geometry, Abort[]];
  writeshp[geometry, filepath];
  writedbf[assoc, StringReplace[filepath, ".shp" -> ".dbf"]];
  filepath
  ]

For point features:

cities = {Entity["City", {"Houston", "Texas", "UnitedStates"}], 
   Entity["City", {"SanAntonio", "Texas", "UnitedStates"}], 
   Entity["City", {"Dallas", "Texas", "UnitedStates"}]};
geometry = 
  Point /@ Reverse /@ (LatitudeLongitude[#] & /@ cities // QuantityMagnitude);
data = <|"Name" -> {"Houston", "San Antonio", "Dallas"}, 
  "Population" -> (CityData[#, "Population"] & /@ cities // QuantityMagnitude), 
  "Elevation" -> (CityData[#, "Elevation"] & /@ cities // QuantityMagnitude)|>

exportshapefile["filepath\\lines.shp", geometry, data];

Now let's pull it in to ArcMap. Since the .prj file isn't a standard format, we don't write it, but you can use the Define Projection toolbox with ArcGIS to do that. It will create the .prj file for us:

enter image description here

enter image description here

We can then just copy the newly created prj file for use in the other examples.

Polylines:

geometry = 
  Line /@ Subsets[Reverse /@ (LatitudeLongitude[#] & /@ cities //  QuantityMagnitude), {2}];
distances = (GeoDistance @@@ Map[Reverse, geometry[[All, 1]], {2}] // QuantityMagnitude);
cesnaspeed = 
  WolframAlpha["speed of a cesna in mph", {{"Result", 1}, "ComputableData"}] // QuantityMagnitude;
data = <|"ID" -> {1, 2, 3}, "Distance" -> distances, "CesnaTime" -> distances/cesnaspeed|>

exportshapefile["filepath\\lines.shp", geometry, data];

enter image description here

And lastly, polygons:

geometry = {Polygon[Reverse /@ (LatitudeLongitude[#] & /@ cities // QuantityMagnitude)]};
data = <|"Name" -> {"The Triangle!"}|>;

exportshapefile["filepath\\poly.shp", geometry, data];

enter image description here

The ESRI format supports other shape types. And it would be trivial to extend the above functions to handle those, I just didn't have a need at the time. Hope this helps someone in the future. Thanks to @WReach's answer here for providing a lot of the knowledge for binary files.

$\endgroup$
  • $\begingroup$ Could you make a package out of this and host it somewhere else than SE? (GitHub?) It looks quite useful. $\endgroup$ – Szabolcs Aug 6 '15 at 21:53
  • $\begingroup$ @Szabolcs I could.... if I knew how to operate GitHub... which I don't--embarrassingly. $\endgroup$ – kale Aug 6 '15 at 21:54
  • $\begingroup$ @Szabolcs Okay, I added this to GitHub. Can you check and make sure I did it correctly. $\endgroup$ – kale Aug 7 '15 at 0:54
  • $\begingroup$ You might want to look at this to "teach" Export[]. $\endgroup$ – J. M. will be back soon Aug 7 '15 at 0:55
  • $\begingroup$ @Guesswhoitis. Interesting. Never knew you could do that. I'll add this later... Right now, if I think anymore about shapefiles, I may stroke out. $\endgroup$ – kale Aug 7 '15 at 0:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.