Jun 17, 2009

Mathematica 7: Export ESRI shapefile

One way to achieve this goal is by using the spatial database. Oracle and PostgreSQL have excellent supports of spatial data types. All you need to do is to connect the database in Mathematica and insert the data into the spatial database. There are tools come with databases allow you to dump the data into a shapefile. If you need a light-weighted spatial database, you can try SpatiaLite, it is based on popular SQLite.

For Windows platform, download the followings first:

spatialite-tools and init_spatialite-2.3.sql

upzip them and copy spatialite.exe and init_spatialite-2.3.sql into the same folder.

In this example, we like to export the following information to a shape file:

data={First[#], CityData[#,"Population"], CityData[#,"Longitude"], CityData[#,"Latitude"]} &/@ CityData [{All, "Indiana", "UnitedStates"}];

We need write some SQLs:

CREATE TABLE Towns (Name TEXT, Population INTEGER);
SELECT AddGeometryColumn('Towns','LonLat', 4326,'Point',2);

4326 means  EPSG 4326, the coordinate reference system of WGS84(longitude, latitude) pair coordinates in degrees.

To Insert the data:

INSERT INTO Towns (Name, Population, LonLat) Values ('Indianapolis', 784118, GeomFromText ('Point(-86.1477 39.7909)', 4326));

In Mathematica, this will create all the “INSERT” statements

str="INSERT INTO Towns (Name, Population, LonLat) Values ('v1', v2, GeomFromText('Point(v3 v4)',4326));"

strs=StringReplace[str, {"v1"->#[[1]], "v2"->ToString[#[[2]]],"v3"->ToString[#[[3]]], "v4"->ToString[#[[4]]]}] &/@ data;

Export["insert.txt", strs]

Then we can put all the SQL statements together into one file (test.sql):

CREATE TABLE Towns (Name TEXT, Population INTEGER);
SELECT AddGeometryColumn('Towns','LonLat', 4326,'Point',2);
… INSERT INTO commands is here …

In window command line mode:

spatialite.exe -init init_spatialite-2.3.sql test.sqlite < test.sql

So far, we have got all the data from Mathematica into the table Towns in test.sqlite.

run: spatialite.exe test.sqlite, and check if the records inside are right, then dump the shapefile:

.dump Towns LonLat towns_shp ASCII POINT;

Import the town_shp.shp back into Mathametica:

Show[Import["towns_shp.shp"], Frame -> True, FrameTicks -> Automatic]


More information: SpatiaLite Tutorial

If you are not comfortable with command-line tools, there is a GUI tool for SpatiaLite.

No comments: