How to make arrow objects for Google Earth

In the previous post, we has use GeoDestination function to make circle objects, and we can use the same function to generate arrow objects, too.

The data we have is the position and displacements in north and east direction, {lat, lon, ndis, edis}. ndis, edis is in millimeter in this case.

First, let calculator the length of the arrow, the scale here is 1 mm displacement = 1000 meters

length = Sqrt[ndis^2 + edis^2]*1000

Then the angle. GeoDestionation requires GeoDirection, it starts from the north.

angle = ArcTan[edis, ndis] /Degree
geoangle = 90 - angle

For the end point:

end = First@GeoDestination[{lat, lon}, {length, geoangle}]

The trick to generate the arrow is that to form another smaller circle around the end point, pick up two points from the arrow heads. 

arr1 = First@GeoDestination[end, {dist*.25, geoangle + 180 - 30}]
arr2 = First@GeoDestination[end, {dist*.25, geoangle - 180 + 30}]

Then we can have two lines which draw in this order, {org –> end}, {arr1->end->arr2}. When comes to export kml, we need to use <MultiGeometry> object in KML, it seems not supported by Export function in Mathematica yet, so you can just export the following string directly:

<MultiGeometry>
<LineString>
<tessellate>1</tessellate>
<coordinates>
-117.093195026,34.116408002,0 -117.020521006,34.0741018797,0
</coordinates>
</LineString>
<LineString>
<tessellate>1</tessellate>
<coordinates>
-117.042636322,34.0757361767,0 -117.020521006,34.0741018797,0 -117.029869517, 34.0907835742, 0
</coordinates>
</LineString>
</MultiGeometry>

Here is the screenshot of the final product in Google Earth:

How to make a circle for Google Earth

KML doesn’t have the circle object built in, a circle can be made by line or polygon object. With GeoDestination function, a perfect circle can be created.

Polygon@Table[GeoDestination[GeoPosition[{39.17, -86.52 }], {5000, a}], {a, 0, 360, 10}][[All, 1, {2, 1}]]

This line will generate a 36-sides polygon centered at 39.17, –86.52 with 5000 meters radius. Then export it to kml, you probably has to manually modify the kml file to set up the color styles.

Here is an example, we use the radius to represent the vertical motion, 1 cm motion = 5000 meters on the ground.

Fun with Texture

Texture support is one major improvement in Mathematica 8.0. Sure, we can have some fun with it. Here is an example. We draw some polygons first.

Then some textures are applied, we get the following image. The main trick is to use alpha channel to get this cut-out texture mapping effect. The original tree image has the alpha channel already, we can use it directly, Texture[ImageData[tree]] will do the trick, if you use Texture[tree], somehow it will not pick up the alpha channel. For the character image, it doesn’t has the alpha channel, we can set a alpha channel by the following command:

SetAlphaChannel[im,
ChanVeseBinarize[ColorNegate@im] // FillingTransform ]

By the way, the character is Tom Nook: "This isn't fat, I've just lined my insides with money."

Download the notebook.

Top 10 Box Office Movies Wallpaper Generator

I saw a script IMPloader that fetches the posters for the top 10 Box Office movies of the week from http://www.impawards.com/, and automatically generates a wallpaper, it is requires Ruby, imagemagick and wget. If you are running Window, you probably don’t want to install these packages. We can do it with Mathematica with ImageAssemble and ImagePad functions introduced in 8.0.

Download the notebook: movieposters.nb

Control Google Earth from Manipulate

One reader asks: could one drive Google Earth from a set of controls within a Manipulate? The answer is yes, however, I only know how to do it with Windows platform through Google Earth COM API. In the past, this issue is covered in this post: Control Google Earth from Mathematica.

The example here is very simple, rotate the camera from Manipulate.

Needs["NETLink`"]
InstallNET[];
(* start up goolge earth *)
ge = CreateCOMObject["GoogleEarth.ApplicationGE"];
(* get the camera object *)
cam = ge@GetCamera[0]
(* define a function to rotate the camera *)
runcam[{lat_, lon_}] := Module[{},
   cam@FocusPointLongitude = lon;
   cam@FocusPointLatitude = lat;
   ge@SetCamera[cam, 2]];

Manipulate[
runcam[{lat, lon}]; {lat, lon}, {lat, -50, 50}, {lon, -180, 180}]

Visualize underground structures with texture mapping

The Center for Remote Sensing of Ice Sheets (CReSIS) has many radar images which show underground ice bed elevations both in north pole and south pole regions.

In the following image, the red line marks the ice bed underground (6000 actually means 6000 m under the earth surface). Radar signal is acquired from the plane, each image has the flight path information (red line in the second image).

The task here is to combine flight path and radar image to display 3D underground information directly. It can be achieved in the following two steps with Mathematica 8’s texture support:

f = BSplineFunction[flightpath]; (* flight path {{x1,y1}, ..{xn, yn}} *)

ParametricPlot3D[Flatten@{f[x], z}, {x, 0, 1}, {z, -6000, 1000},
BoxRatios -> 1, BoundaryStyle -> Green, Mesh -> None, Axes -> None,
PlotStyle -> {Texture[radarimage]}, Lighting -> "Neutral", ImageSize -> 600]

The radar image has to be cropped in advance, and we only keep part as deep as 6000 meters. There are some rough parts on the image, especially close to the boundaries (check out the update for the solution), the general quality is good enough. We can make the same thing for another radar image, then put them together, make the following product to show the intersection. It is not perfect solution, however, it is nice to have something in just several lines. 

No notebook this time.

Updates: Someone asks how to get rid of zig-zag boundaries.

It can be done with two options, one is PlotPoints, the other MaxRecursion, check the ParameterPlot3D documentation for the details.

For example,  MaxRecursion –> 4

Pi Day

Reading: A history of Pi

I found the old notebook to generate Pi T-shirt design.

Download Notebook

Mathematica and Spatial Database

Warning: this post uses undocumented Mathematica command and modifies the Mathematica installation. Don't try this at home.

For any GIS software, evaluating spatial relationships, such as equal, disjoint, within, intersects, etc., is a fundamental requirement. Also, R-tree support is a must for any large spatial data set. Currently, these sets of functions are not built in Mathematica. It is difficult to perform complex GIS analysis inside Mathematica. One way is to call external libraries through Mathlink/Jlink. Another way is connecting to a spatial database, such as PostGreSQL, Oracle Spatial through database connection. Spatial database? What about Spatialite, it is a complete spatail DBMS built as an extension to the extreme light-weighted database SQLite. 

From this Wolfram|Alpha tweet analysis post, it shows Mathematica actually ship with SQlite.

I installed the Mathematica 8 linux trail version.

db = Database`OpenDatabase["/tmp/test.sqlite"];
Database`QueryDatabase[db, "SELECT sqlite_version();"]
{{3.6.1}}

Here is SQLite library from Mathematica
/usr/local/Wolfram/Mathematica/8.0/SystemFiles/Kernel/Binaries/Linux/libsqlite3.so

For security reason, the dynamic loading extension is disabled. And we need to compile our own copy of libsqlite3.so. The detail is explained here: the pre-packaged 'libsqlite' trap.

Grab the source code for SQLite website, build the new library with:

CFLAGS="-DSQLITE_ENABLE_LOAD_EXTENSION=1" ./configure

then make a copy of the original libsqlite3.so, and overwrite it with the new version.

Database`QueryDatabase[db, "SELECT sqlite_version();"]
{{3.7.5}}

Then we can try with the Spatialite already installed:

Database`QueryDatabase[db, "SELECT load_extension(‘/usr/lib/libspatialite.so’);"]
Database`QueryDatabase[db, "SELECT spatialite_version();"]
{{"2.4.0"}}

Try a spatial SQL command:

Database`QueryDatabase[db, "SELECT X(GeomFromText('POINT(-85 39)',4326));"]
{{-85.}}

Wow, it is working!

Here is a tutorial on Spatialite, you can get the idea of what kind of functions are supported by Spatialite. If you are familiar with the spatial database, you can build a rather functional GIS system inside Mathematica with the support of Spatialite.

Create simple DEM from Google Map

With texture support in Mathematica 8.0, we can create the simple DEM with image overlay by combing the functions from two previous posts Google Static Map and Google Elevation API.

 

However, there is the usage limits on Google Elevation API. You’d better pull the dem once and export the data, otherwise, you may hit the daily limits very quickly.

Here is an example on displaying earthquake data with the dem, the red line outlines the fault plane.

Download notebook.