Jun 19, 2009

Control Google Earth from Mathematica

This is for Windows platform only. With Mathematica’s .NET link and Google Earth COM API, we can control Google Earth’s camera and add features directly from Matheamtica.

There is an example of planning a shortest tour through every country of the world in Document: FindShortestTour.

SC[{lat_,lon_}]:=r {Cos[lon \[Degree]] Cos[lat \[Degree]],Sin[lon \[Degree]] Cos[lat \[Degree]],Sin[lat \[Degree]]};
distfun[{lat1_,lon1_},{lat2_,lon2_}]:=VectorAngle[SC[{lat1,lon1}],SC[{lat2,lon2}]] r;

Let’s view this example in Google Earth:


(* startup google earth *)
ge = CreateCOMObject["GoogleEarth.ApplicationGE"];

(* load path file already generated *)

(* get the camera object *)

(* funcion to ratate the camera *)

(* let's see the movie *)
runcam[#]&/@ centers[[route]];

Here is a low quality video.


Here is test2.kml for the shortest path

Here is the recorded tour (shortestrout.kmz) in Google Earth

First load test2.kml to Google Earth, then double click ShortestRoute.kmz to view the animation.

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.