Can Mathematica Export Shapefiles?

Question

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

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

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

Answer 1

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, {5, 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:

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

And lastly, polygons:

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

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

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.

Answer 2

Thanks for making a very useful Package even in 2021. However, Polygon output seems to have an error, checked with QGIS 2.18. I found the error part in writepolygon[ ] function for entering Shape Type value(Polygon:5) as follows.

[Original code]

writepolygon[stream_, shxstream_, polyrecord_, recordnumber_] := 
 Module[{numpart = Depth@polyrecord[[1]] - 2, numpoints, bounds},
 .... 
  BinaryWrite[
   stream, {3, bounds[[1, 1]], bounds[[2, 1]], bounds[[1, 2]], 
 ...

[Revised code]

writepolygon[stream_, shxstream_, polyrecord_, recordnumber_] := 
 Module[{numpart = Depth@polyrecord[[1]] - 2, numpoints, bounds},
 .... 
  BinaryWrite[
   stream, {5, bounds[[1, 1]], bounds[[2, 1]], bounds[[1, 2]], 
 ...