Lunchtime Playground: Fun with Mathematica

How to make arrow objects for Google Earth

2011-10-31T11:29:00.001-04:00

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:

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

How to make a circle for Google Earth

2011-06-02T10:00:00.001-04:00

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

2011-05-11T15:00:00.001-04:00

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

2011-05-09T08:51:00.001-04:00

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

2011-04-29T12:31:00.001-04:00

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

2011-04-26T10:20:00.001-04:00

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

2011-03-14T09:33:00.001-04:00

Reading: A history of Pi

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

Download Notebook

Mathematica and Spatial Database

2011-03-11T12:45:00.001-05:00

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

2011-02-04T13:15:00.001-05:00

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.

Happy Holiday!

2010-12-21T09:55:00.001-05:00

Happy Holiday!

Be a nice person in 2011, everyone will be happy.

How to generate this image?

First, find an image, then

img2 = ColorConvert[img, "HSB"];
img3 = ImageData[img2][[All, All, 3]];
ListPlot[Table[img3[[i]]*10 - i, {i, 1, ImageDimensions[img][[2]], 3}],
PlotRange -> All, Axes -> None, Joined -> True, PlotStyle -> {{Black}, {Thick, Lighter[Black]}}, AspectRatio -> ImageAspectRatio[img]]

Color-coded contour lines with Mathematica

2010-11-15T11:05:00.001-05:00

In ContourPlot, ContourStyle –> function can be used to color-code contour lines.For example, the number of contours is numC=20, then for each contour, the color is defined by a coloring function:

ContourStyle -> Table[{ColorData["Rainbow", (i - 1)/(numC - 1)]}, {i, numC}]

Of course, it is also necessary to set Contours –> numC.

Here are the examples, the famous peaks function from Matlab is used.

2D ContourPlot example:

Turn 2D plot into the 3D one:

Another 3D example:

Download the notebook.

Customizing DateListPlot with PlotMarkers

2010-11-08T10:07:00.001-05:00

This example shows you how to specific PlotMarkers for each point in a data set. The solution is quite simple, just partition the dataset into different series, each set only contains one point exactly: Partition[data, 1]. Then we can assign the different PlotMarkers for each point.

Let’s download some weather data:

data = WeatherData[$GeoLocation, "WindDirection", {{2009, 1, 1}, {2009, 1, 5}}];
data = Select[data, FreeQ[#, {_, Missing["NotAvailable"]}] &];

We like to use arrows to represent wind directions:

markers =
Graphics[{Red, Arrow[{{0, 0}, -0.5 {Sin[#[[2]] Degree], Cos[#[[2]] Degree]}}],
EdgeForm[], FaceForm[], Rectangle[{-1, -1}, {1, 1}]}] & /@ data;

We put a invisible box around arrow to make sure that the markers are aligned by {0, 0}.

newdata = data; newdata[[All, 2]] = 0;
g = DateListPlot[Partition[newdata, 1], PlotRange -> All, PlotMarkers -> markers, Axes -> {True, False}, FrameTicks -> {Automatic, None}]

Let’s test another place: Tokyo, Japan

Summer time:

Winter time:

The pattern of winter “north” wind and the summer “south” wind is very clear, it is much better than just plotting points.

No notebook, all the codes are here already.

Plotting on Google static map

2010-11-01T12:53:00.001-04:00

We have given an example on plotting with WMS. You may be wondering how we can do the similar thing with Google static map. It is very easy to use with Mathematica:

googleMap[{lat_,lon_},{sizex_,sizey_},zoom_,maptype_]:=Import["http://maps.google.com/maps/api/staticmap?sensor=false&center="<>ToString[lat]<>","<>ToString[lon]<>"&zoom="<>ToString[zoom]<>"&size="<>ToString[sizex]<>"x"<>ToString[sizey]<>"&maptype="<>maptype];

Google map is based on Mercator Projection. So the data has to be projected, rather than drawing (longitude, latitude) pairs directly. The scale on east-west (longitude) is constant with zoom level, and only latitude needs to be projected.

(* convert lat to y *)
lat2y[lat_]:=0.5Log[(1+Sin[lat Degree])/(1-Sin[lat Degree])];
(* convert y to lat *)
y2lat[y_]:=Gudermannian[y]*180/Pi;

We have the data in pairs of (longitude, latitude) which lie in a bounding box ((lon1,lat1), (lon2, lat2)), and to plot them on a Google static map, we need: 1. find the center of the map; 2. find the right zoom level; 3. find the proper map size; 4. project the data: (longitude, lat2y[latitude])

I will skip the details, if you are interested into the mathematics, check out these two posts, : R-Google Map and How to make Google static maps interactive.

Here are two examples:

coords = CountryData["Australia", "Coordinates"];
(* only project latitude *)
Map[{#[[2]],lat2y[#[[1]]]}&, coords,{2}]]

In the second example, GIS data is imported from an XML file, it is from the post before Mathematica 7.0 released.

Mathematica 8 has the texture function, so we can import static Google map as the texture, it will be easy to create complex 3D GIS visualization.

Download the notebook for detail.

Processing comic images for e-reader

2010-10-26T17:20:00.001-04:00

First, this post probably falls into the category of “just because you can do it doesn’t mean you should” with Mathematica.

Here is the job: the comic is in cbz/cbr format, and images are too large for the e-reader usually with 800 by 600 resolutions (the actual usable space is even smaller). So you like to do three things with these images: 1.cut white margins, 2. resize, 3. increase contrast if image is not clear enough. It can be done with three Mathematica commands: ImageCrop, ImageResize and ImageAdjust.

One line example:

Export[#,
ImageAdjust[
    ImageResize[
    ImageCrop[Import["test.cbz", {"ZIP", #}]], {584, 714}],
    0.2]] & /@ (Import["test.cbz", "ZIP"])

If it is in cbr format, it has to be unpacked first.

Actually under the certain situation, this may come in handy. There are more and more comics in epub format made for iPad. And epub is basically a zip file contains html/css and images. The following line returns the image file names inside epub.

Select[Import["test.epub","ZIP"], StringMatchQ[#,__~~".jpg", IgnoreCase->True]&]

Then the images can be changed into grayscale if necessary and re-formatted for e-reader. Just dump the images back into epub afterward.

Tips: Import MATLAB MAT-files

2010-09-29T15:58:00.001-04:00

You want to import the “.mat” file:

Import[“Helheim.mat”]

Import::fmterr: Cannot import data as MAT format. >>

Don’t panic. Since the version 7.3 (2006b), Matlab actually saves mat files using the HDF5 format by default. So pretty much every mat file you work with nowadays is in HDF5 format. Let’s try it again.

Import["Helheim.mat", "HDF5"]

{"/Data", "/Depth", "/Distance", "/Latitude", "/Longitude", "/UTC_time"}

data = Import["Helheim.mat", {"HDF5", "Datasets", "/Data"}];

Then we can generate a beautiful plot in Mathematica.

For who may wonder what this image is, it is the radar image shows underground ice-thickness, more information can be found at The Center for Remote Sensing of Ice Sheets (CReSIS)

Work with Web Map Service (WMS)

2010-08-25T20:47:00.001-04:00

Web Map Service (WMS) is is a standard protocol for serving geo-referenced map images over the Internet. JPL OneEarth is used as the data source here.

All the datasets that exist on OneEarth server are documented for WMS clients in an XML Capabilities file, which contains information about the layers, which correspond to datasets, and the styles associated with them.

The base url is http://wms.jpl.nasa.gov/wms.cgi?, to request the image, some parameters are required: layers, srs, width, height, bbox, format, styles. Click the following url, it returns an satellite image from Blue Marble Next Generation, a MODIS-derived 500m true color earth dataset showing seasonal dynamics. styles is used to specify the month. Check OneEarth website for the details.

http://wms.jpl.nasa.gov/wms.cgi?REQUEST=GetMap&layers=BMNG&srs=EPSG:4326&width=300&height=300&bbox=112,-44,154,-10&format=image/jpeg&styles=Jan

In Mathematica, image is imported with Import[wms_url]. However, once the image is imported, it loses the geo-reference information, and we need move the image back to its right position (bbox parameter).

The trick is using Raster function to convert image to a graphic object with the right coordinates.

raster = Raster[Reverse[ImageData[img]], {{112, -44}, {154, -10}}];

I borrow the example from this weatherdata blog entry.

The imported image works perfectly with the data from Mathematica.

Mathematica Notebook: wmsimage.nb

Visualize weather pattern with random walk

2010-08-20T12:44:00.001-04:00

When talking about the “good” weather, the daily temperature change is important, slow steady change from day to day is much better than sudden temperature increase or drop.

Here is the daily mean temperature for Bloomington, Indiana.

We define a piecewise function to classify the day to day temperature change. x is (T – Tprevious day), y is the threshold to define the weather change: 1 ~ no change, 2 ~ cold move, 3 ~ mild cold, 4 ~ mild hot, 5 ~ hot move.

Then we can apply the same random walk trick used in the post Visualize irrational number as random walk. The classification results are mapped into the moves: 1 ~ no move {0, 0}, 2 ~ move south {0, –1}, 2 ~ move west {-1, 0}, 3 ~ move east {1, 0}, 4 ~ move north {0, 1}. Then we can get a random walk image (the green dot is starting point {0, 0}, the red one is the end).

This is probably nothing interesting. However, we can use it to compare the patterns among different cities. Here is the results from Washington, Columbus, Indianapolis, Kansas City, Denver. The threshold is 3.5 degree. This shows the certain pattern among cities from ocean to inland.

When we compare the cities across the ocean, a different pattern is shown.

Down the weatherwalk.nb for the detail.

Play with Bibliographic Data

2010-04-19T23:19:00.001-04:00

I’ve got the publication list of a research institute. There are more than 1400 entries, exported in xml format from Endnote. As a data mining project, there are plenty of things you can do with the data. Let’s assume we are interested in the relationships among researchers and groups. Can we check this out quickly in Mathematica?

First, extract the data:

in the xml, for each publication, the author list is stored as the following:

xml = Import["Bib.xml"];

authors = Cases[xml, XMLElement["authors", _, authors_] -> authors, Infinity];

names = Flatten@Cases[#, XMLElement["author", _, {___, XMLElement["style", _, name_]}] –> name] & /@ authors;

What we get here is the lists of authors for each publication.

Let’ see the relationship between number of authors and number of publications.

Sort[Tally[Length[#] & /@ names], #1[[1]] < #2[[1]] &]

{{1, 265}, {2, 280}, {3, 320}, {4, 224}, {5, 94}, {6, 85}, {7,
39}, {8, 22}, {9, 20}, {10, 14}, {11, 11}, {12, 14}, {13, 5}, {14,
3}, {15, 5}, {18, 1}, {19, 1}, {20, 1}}

We can see that most of publications have no more than 4 authors.

Next step, we like to check out the internal relationships among authors. We need to generate a network for authors. For example, if a publication has 4 authors {A, B, C, D}, the network is defined as a circle:

Flatten@(Partition[Append[authors, authors[[1]]], 2, 1] /. {x_, y_} :> {x -> y})

{A -> B, B -> C, C -> D, D -> A}

For all the publications, we get the following network:

There are two large research groups inside this institute. Then re-draw the graph with top 10 contributors, it confirms the information.

No too bad with 10 minutes coding.

I will not release the notebook this time, since I may not have the right to distribute the data. Sorry about it.

Extract elevation data with Google Elevation Service

2010-03-25T15:03:00.001-04:00

In the previous post: Extract elevation data from Google Earth, Google Earth COM API is used, it only works on Windows platform. Google Map now has Elevation Web Service, it is quite easy to do it with new API. The new service does not require a Maps API key. The basic form is

“http://maps.google.com/maps/api/elevation/outputFormat?parameters”

For output format,

/json returns results in JavaScript Object Notation (JSON).
/xml returns results in XML, wrapped within a <ElevationResponse> node.

JSON is easy to parse, this is a sample query:
json=Import["http://maps.google.com/maps/api/elevation/json?locations=39.7391536,-104.9847034&sensor=false"]

{
  "status": "OK",
  "results": [ {
    "location": {
      "lat": 39.7391536,
      "lng": -104.9847034
    },
    "elevation": 1608.8402100
  } ]
}

ToExpression@StringCases[json, NumberString]

{39.7392,-104.985,1608.84}

Here is the example output:

Path elevation example:

DEM + Path:

Grab the GoogleElevationService.nb for detail.

Update: answer the comments

You may have the trouble with the notebook. Sometimes ListPlot3D runs too slow and even crashes Mathematica. You can switch to ListPlotPoint3D.

ListPointPlot3D[dem[[All, {2, 1, 3}]], ColorFunction -> "Rainbow"]

Here is the example with 100 by 100 DEM.

More on Heatmap

2010-01-26T19:16:00.001-05:00

A reader asks about this heatmap post on Flowing Data. Sure, we can do it easily with ArrayPlot. Grab the notebook here.

data = Import["ppg2008.csv"];

Grid[data]

numbers = data[[2 ;; All, 2 ;; All]]; (* this is the data for heatmap *)
playernames = data[[2 ;; All, 1]]; (* for the labels *)
statnames = data[[1, 2 ;; All]]; (* for the labels *)

Then we need to scale the each column separately to [0,1].

newnumber = Transpose[Rescale[#] & /@ (Transpose[numbers])];

Then you can generate a simple heatmap simply by ArrayPlot[newnumber].

Of course, we like to add the labels with FrameTicks option. FrameTicks->{{left,right},{bottom,top}} mark options specified separately for each edge.

Transpose[{Range[Length[playernames]], playernames}] generate the list {{1, “player1”} … {50, “player50”}} to label the player names.

Transpose[{Range[Length[statnames]], Map[Rotate[#, 60 Degree] &, statnames]}] with this line, we can rotate the labels at the same time.

Here is the final product, click to see the full graphic.

Charting time series as calendar heat maps

2010-01-12T12:26:00.001-05:00

I haven’t updated this blog for a while. I am working on one Matlab project right now. It is kind of difficult for me to work with Matlab and Mathematica at the same time.

There is an interesting post Charting time series as calendar heat maps in R. The original idea from SAS Analysis of airline performance. I create a simple version in Mathematica. The tricky part is to generate the background grid.

The code I use is from an old tutorial of Making a Calendar. For each month, the boundary is defined by 8 points, marked out by darker line. The monthly grids are shifted to the right positions to form a yearly grid. Once you figure out this part, the rest is just straightforward.

I use the stock as the example, this is actually not the best data set for this type of visualization.

Calendar heatmap:

Download the file calendarheatmap.nb for the detail.

By the way, you can do whatever you like with the materials posted on this blog, there is no copyright problem.

Make a US county thematic map using Mathematica

2009-11-13T16:38:00.001-05:00

I read a blog entry How to Make a US County Thematic Map Using Free Tools. The idea is to use Python script to modify the style settings in SVG file to make the most recent unemployment map. Mathematica can’t import SVG file directly, however, SVG file is in XML format, it is very easy to extract the necessary data from SVG file. If you are not sure what’s going here, please refer to the original blog.

All the files are zipped together, just unzip it, run the notebook. Download it here.

User-defined color themes

2009-11-02T16:02:00.001-05:00

With Blend function, it is quite simple to use user-defined color themes.

Blend[{col₁, col₂, col₃, ...}, x]: linearly interpolates between colors col_i as x varies from 0 to 1.

We like to use the following color theme:

c = {{37, 57, 175}, {40, 127, 251}, {50, 190, 255}, {106, 235,
255}, {138, 236, 174}, {205, 255, 162}, {240, 236, 121}, {255,
189, 87}, {255, 161, 68}, {255, 186, 133}, {255, 255, 255}};

colors = RGBColor[#/255] & /@ c;

This shows the each color in the theme:

Graphics[Table[{EdgeForm[Black], FaceForm[colors[[i]]],
Rectangle[{i, 0}, {i + 1, 1}]}, {i, 1, Length[colors]}]]

Check the color theme with Blend function:

DensityPlot[x, {x, 0, Length[c]}, {y, 0, 1}, AspectRatio -> Automatic,
FrameTicks -> None, ColorFunction -> (Blend[colors, #] &),
PlotRangePadding -> None]

Then you probably notice how to use it in your own plot,

ColorFunction -> (Blend[colors, #] &)

Test the color theme with the data:

ReliefPlot[data, ColorFunction -> (Blend[colors, #] &)]

Maybe the ligher color is better.

lightercolors = Lighter[#] & /@ colors;

Just for fun, let’s play the color theme with an existing image.

img=ImageData[ColorConvert[place_any_image_here, ”Grayscale”]];

ArrayPlot[img, ColorFunction -> (Blend[darkercolors, #] &)]

Wikipedia Page Analysis

2009-10-20T13:20:00.001-04:00

Wikipedia has lots of scientific information, however, due to its nature, it is still not considered as a research resource. This doesn’t mean it has to be ignored. I have checked some pages related with various topics in GIS field. Most of them are well-written, the information are actually quite accurate, several contributors are the professionals in the field. In this post, I like to check some metadata information of “Mathematica” Page on Wikipedia, it may gives us some ideas about its quality.

Tools we need: Mediawiki API and Mathematica. There are plenty examples on how to use Mediawiki api. Basic procedure is to use Import[queryurl,”XML”], then parse xml to get the information we need.

Page revision history:

(* import contributor and timestamp *)

url = "http://en.wikipedia.org/w/api.php?action=query&prop=revisions&\
titles=Mathematica&rvprop=user|timestamp&rvlimit=500&redirects$rvuser&\
format=xml";

xml = Import[url, "XML"];
rawdata= Cases[xml, XMLElement["rev", w_, _] :> w, Infinity];
data = {"user", "timestamp"} /. rawdata;

This page is constantly revised, we probably can assume the information on “Mathematica” page is up-to-date.

The information on the contributors is also interesting.

We can dig out more information on the contributors:

(* import paged edited by each user *)

userpages[usr_] :=
Module[{url, uxml, udata, unicase},
   url = "http://en.wikipedia.org/w/api.php?action=query&list=\
usercontribs&uclimit=500&format=xml&ucuser=" <> usr;
   uxml = Import[url, "XML"];
   udata = Cases[uxml, XMLElement["item", w_, _] :> w, Infinity];
   unicase = DeleteCases[Union["title" /. udata ],
     x_ /; (StringMatchQ[x, "User talk:" ~~ __] || StringMatchQ[x, "Talk:" ~~ __] || StringMatchQ[x, "User:" ~~ __])]; Map[usr -> # &, unicase]];

The common pages edited by these top5 contributors:

From the pages they have edited, they have worked on several topics closely related with Mathematica. This looks good, we may say they probably know what they are doing.

Update:

Download Wikipedia Notebook for the details.

Mathematica 7 on Windows 7

2009-10-15T23:46:00.001-04:00

We just updated our machine from Windows Vista to Windows 7. I re-installed Mathematica 7 and tested with several notebooks. I didn’t notice any difference on performance.

You can re-use the license file from Vista.

C:\Users\user_name\AppData\Roaming\Mathematica\Licensing\mathpass

One feature I tested is Handwritten Math Recognition in Windows 7. It works very well with some simple on-screen drawings. If you use a tablet PC, this may be quite convenient in classroom.

Monitor Mathematica computation with email

2009-09-23T12:19:00.001-04:00

One of the new features in Mathematica 7.0 is allowing the sending of email directly from any Mathematica program. It turns out very convenient for the certain situation. We have some computation tasks take long time to finish, so we let it run on the server. The problem is that once a while I have to ssh to the server to check the outputs to make sure nothing wrong with the computation. Now with the email function, I can get the update immediately for each computation step. This is really helpful, especially now days you can check the email almost everywhere.

If you don’t need to check the detail, sending the update through the twitter is probably more cool.

One thing I am not sure is that how to catch the error. If there is an error raised during the computation, the math kernel can automatically send out the email with the error message, this function will be perfect.

By the way, this week I am working on an algorithm related with the Traveling Salesman Problem. If you check TSP on MathWolrd, download the notebook, you will notice it is created by a newer version of Mathematica. Open the file with any text editor, you will see:

(* CreatedBy='Mathematica 8.0' *).

Visualize irrational number as random walk

2009-08-25T22:26:00.001-04:00

Irrational numbers have decimal expansions that neither terminate nor become periodic. So we can get unlimited “random walk” steps from an irrational numbers.

With the following code, the first 10000 digits of Sqrt[2] is presented as a random walk by converting it in base 4. 0, 1, 2 and 3 digit in base 4 represent 4 directions. Starting point is the green dot, the red one is the ending point.

x = N[Sqrt[2], 10000];
walk = First@ RealDigits[x, 4];
rn = FoldList[Plus, {0, 0}, {{0, 1}, {1, 0}, {0, -1}, {-1, 0}}[[# + 1]] & /@ walk];
Graphics[{Line[rn], PointSize[Large], Green, Point[First@rn], Red, Point[Last@rn]}]

We can also display it with ArrayPlot by constructing a sparse array.

(* shift the moves to {1,1} *)

minx = Min[rn[[All, 1]]];
miny = Min[rn[[All, 2]]];
m = # + {-minx + 1, -miny + 1} & /@ rn;

(* sparse array *)

tt = Tally[m];
cd = SparseArray[tt[[All, 1]] -> tt[[All, 2]]];
cd = Transpose[cd];
ArrayPlot[cd, ColorFunction -> "Rainbow", DataReversed -> True, ColorRules -> {0 -> White}]

Let’s visualize Sqrt[2], e, Pi in their first 50000 digits. It seems there are some similarities among these images.

Sqrt[2]

Update:

DrMajorBob said... in comments:

Here's a 3D version:

x = N[Sqrt[2], 10000];
walk = First@RealDigits[x, 6];
rn = FoldList[ Plus, {0, 0, 0}, {{0, 1, 0}, {1, 0, 0}, {0, -1, 0}, {-1, 0, 0}, {0, 0, 1}, {0, 0, -1}}[[# + 1]] & /@ walk ];
Graphics3D[{Line[rn], PointSize[Large], Green, Point[First@rn], Red, Point[Last@rn]}]

You can also try non-irrational number, e.g. 121/5^10

Thanks.

Update 2:

Visualize genome sequence. I know nothing about it, it probably totally meaningless.

GenePlot[g_] := Module[{cs, rn}, walk = Characters[GenomeData [g, "FullSequence"]] /. {"A" -> 1, "T" -> 2, "G" -> 3, "C" -> 4};
rn = FoldList[Plus, {0, 0}, {{0, 1}, {1, 0}, {0, -1}, {-1, 0}}[[#]] & /@ walk];
Graphics[{Line[rn], PointSize[Large], Green, Point[First@rn], Red, Point[Last@rn]}, Frame -> True, FrameTicks -> None] ]

GenePlot["DDK4"]

GenePlot["DKK3"]

Display a plot by clicking a button

2009-08-12T11:30:00.001-04:00

Ok, this seems to be very easy in Mathematica.

Button["Show me a plot", Show[Plot[Sin[x], {x, -Pi, Pi}]]]

Then you click the button, nothing happens.

Button["Show me a plot", Print[Plot[Sin[x], {x, -Pi, Pi}]]]

This line will do the job. Or use the following line:

Show[Plot[Sin[x], {x, -Pi, Pi}],
DisplayFunction -> (Button["Show me a plot", Print[#]] &)]

The key here is to use “Print” rather than “Show”. It isn’t clearly explained in the “ref/Button”.

View weighted graph with GraphPlot

2009-08-07T15:59:00.001-04:00

Here is a simple example on how to customizing Graphplot. We like to use GraphPlot to visualize the number of people who commute into or out Monroe county from/to its neighbor counties.

g={{"Owen" -> "Monroe", 2813}, {"Greene" -> "Monroe",
3788}, {"Lawrence" -> "Monroe", 4022}, {"Jackson" -> "Monroe",
85}, {"Brown" -> "Monroe", 689}, {"Morgan" -> "Monroe",
821}, {"Monroe" -> "Owen", 676}, {"Monroe" -> "Greene",
207}, {"Monroe" -> "Lawrence", 679}, {"Monroe" -> "Brown",
303}, {"Monroe" -> "Morgan", 617}}

vercoor={"Monroe" -> {-86.529, 39.1621}, "Owen" -> {-86.7642, 39.2868}, "Greene" -> {-86.9403, 39.0246}, "Lawrence" –> { -86.4923, 38.8627}, "Jackson" -> {-86.0462, 38.8798}, "Brown" -> {-86.2382, 39.203}, "Morgan" -> {-86.4238, 39.4233}}

First try:

GraphPlot[g, VertexLabeling -> True, VertexCoordinateRules -> vercoor]

Using arrow to indicate in/out seems to be a good idea. We use EdgeRenderingFunction in second try:

GraphPlot[g, VertexLabeling -> True,
EdgeRenderingFunction -> (Arrow[#1, 0.05] &),
VertexCoordinateRules -> vercoor]

However, the labels on the edge is lost. We can handle it in EdgeRenderingFunction.

GraphPlot[g, VertexLabeling -> True,
EdgeRenderingFunction -> ({Text[#3, Mean[#1]], Arrow[#1, 0.05]} &), VertexCoordinateRules -> vercoor]

The graph is still difficult to read, the commuting pattern isn’t clear at a glance. We further update EdgeRenderingFunction and use the line color and thickness to show the pattern.

GraphPlot[g,
EdgeRenderingFunction -> ({If[#2[[1]] == "Monroe", Red, Blue],
AbsoluteThickness[0.5 + #3/500], Arrowheads[0.02 + #3/120000], Arrow[#1, 0.05]} &), VertexLabeling -> True,
VertexCoordinateRules -> vercoor]

In the last try, we use VertexRenderingFunction to make the label more clear.

GraphPlot[g,
EdgeRenderingFunction -> ({If[#2[[1]] == "Monroe", Red, Blue],
AbsoluteThickness[0.5 + #3/500], Arrowheads[0.02 + #3/120000], Arrow[#1, 0.06]} &), VertexLabeling -> True,
VertexCoordinateRules -> vercoor,
VertexRenderingFunction -> ({Text[Style[#2, 14, Bold], #2 /. vercoor, Background -> White]} &)]

Import the shapefile, then you get a map:

View stock market seasonality

2009-07-29T14:06:00.001-04:00

A reader asks about the seasonal chart. Here is a simpler version.

The seasonality is defined as the percent change from the beginning of year to the end of each month.

monthlychange[stock_, year_] :=
Module[{data, firstdateprice, endofmonthprice},
data = FinancialData[stock, "Close", {{year, 1, 1}, {year, 12, 31}}];
data = Map[Flatten[#] &, data]; firstdateprice = data[[1, 4]];
endofmonthprice = Table[Last@Select[data, #[[2]] == i &][[All, 4]], {i, 12}]; (endofmonthprice/firstdateprice - 1)*100]

This will give you the monthly change of AAPL in the year of 2008.

monthlychange["AAPL", 2008]

For multiple years, the average seasonality:

mc = monthlychange["AAPL", #] & /@ Range[2000, 2008];

Mean@mc

For APPL, the following graph shows the average seasonality from 2000 – 2008. According to this graph, it is probably a good idea to hold the stock through Oct and Nov before selling it.

This is a graph shows the seasonality year by year.

View stock market with heatmap

2009-07-28T15:14:00.001-04:00

There is a heatmap tool which let you view the stock price change at a glance. It can be recreated with Mathematica in several lines of code.

Let’s pull the data by format {“symbol”, price change}

stockdata={#,(FinancialData[#]/FinancialData[#,"Close"]-1)*100 }& /@FinancialData["^DJI","Members"];

Then we need to represent the price change by color: Green means up, Red mean down, and the deeper the color, the bigger the change.

max = Max[Abs[stockdata[[All, 2]]]]; (* max change *)

GraphicsGrid[Partition[Graphics[{If[#[[2]]>=0.0,Blend[{White,Green},Rescale[#[[2]], {0,max}]], Blend[{Red,White}, Rescale[#[[2]],{-max,0}]]], Rectangle[], Black,Text[Style[#[[1]]<>"\n"<>ToString[NumberForm[#[[2]],{3,2}]]<>"%", Medium,Bold], {0.5,0.5}]}] &/@ stockdata, 6 ]]

The key function is if the price change is >0, then rescale the change in range (0, max) and get it’s color in Blend[{White, Green}]; if the price change is <0, then rescale the change in range (-max, 0) and use Blend[{Red, White}] to get the right color.

In the ascending order:

We can try other representations, too. For example, we can use the size of disk to represent the change.

Read binary file with ByteOrdering option

2009-07-16T16:12:00.001-04:00

This is probably a very trivial problem for the ones who work on both PC and Unix platforms. I download a small DEM (478 rows by 399 columns) stored as Int16 binary from The Gloabal Land 1-km Base Elevation Project. The range of elevation is 99 ~ 397 meters.

Then in an Apple Power Mac G5:

data = BinaryReadList["mydata.bin", Table["Integer16", {399}], 478];

ArrayPlot[data, PlotRange -> {99, 397}]

Obviously, something wrong with the numeric value loaded from the binary file. It turns out the binary is using little endian. However, the Power Mac is using big endian. You can check your system with $ByteOrdering.

Read the data again with the right ByteOrdering option.

data = BinaryReadList["mydata.bin", Table["Integer16", {399}], 478, ByteOrdering->-1];

ArrayPlot[data, ColorFunction->"Rainbow", PlotRange -> {99, 397}]

Extract elevation data from Google Earth

2009-07-09T10:51:00.001-04:00

In Google Earth COM API, there is a function: GetPointOnTerrainFromScreenCoords

Given an screen_x and screen_y, it returns IPointOnTerrainGE object, which gives out the {Latitude, Longitude, Altitude}

Screen coordinates range from (-1, –1) to (+1, +1)

(-1, -1) - bottom left hand corner of the screen. (0,0) - center of the screen. (1, 1) - top right hand corner of the screen.

Let’s use this function in Mathematica to extract elevation data.

First, zoom Google Earth into a testing place

Then in Mathematica:

Needs["NETLink`"]
InstallNet[]
ge = CreateCOMObject["GoogleEarth.ApplicationGE"];

getAltitude[x_, y_] :=
Module[{pv},
pv = ge@GetPointOnTerrainFromScreenCoords[x, y]; {pv@Longitude, pv@Latitude, pv@Altitude}];

(* extract 50 by 50 grids around the center of the screen *)
dem = Table[
getAltitude[x, y], {x, -0.5, 0.5, 0.02}, {y, -0.5, 0.5, 0.02}];

Next situation: Given a list of {Latitude, Longitude}, how can we get the corresponding elevations for each location?

The tip is to use the camera control to move the center of screen to the given {Latitude, Longitude}.

cam = ge@GetCamera[1];

getAltitudebyLatLon[{lat_, Lon_}] :=
Module[{pv}, cam@FocusPointLongitude = Lon;
cam@FocusPointLatitude = lat; cam@Range = 8000; ge@SetCamera[cam, 5];
pv = ge@GetPointOnTerrainFromScreenCoords[0, 0]; pv@Altitude]

cam@Range is the control of “eye alt”, it is in meters.

It comes very handy for extracting the cross-profile.

Control Google Earth from Mathematica

2009-06-19T17:50:00.001-04:00

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]]};
r=6378.7;
places=CountryData["Countries"];
centers=Map[CountryData[#,"CenterCoordinates"]&,places];
distfun[{lat1_,lon1_},{lat2_,lon2_}]:=VectorAngle[SC[{lat1,lon1}],SC[{lat2,lon2}]] r;
{dist,route}=FindShortestTour[centers,DistanceFunction->distfun];

Let’s view this example in Google Earth:

Needs["NETLink`"]
InstallNET[];

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

(* load path file already generated *)
ge@OpenKmlFile["d:/temp/test2.kml",1]

(* get the camera object *)
cam=ge@GetCamera[1];

(* funcion to ratate the camera *)
runcam[{lat_,lon_}]:=Module[{},
cam@FocusPointLongitude=lon;
cam@FocusPointLatitude=lat;
ge@SetCamera[cam,2]];

(* 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.

Mathematica 7: Export ESRI shapefile

2009-06-17T15:48:00.001-04:00

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):

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

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.

Testing the new map with Mathematica

2009-03-24T12:45:00.001-04:00

Generally speaking, when dealing with GIS data, you probably shouldn’t put Mathematica in top of your list. However, when you want to develop a new algorithm or visualization, Mathamtica is a good choice to test the prototype.

Here is an example:

The same map in the “forest” looking.

Color Blend

2009-03-09T20:49:00.001-04:00

Rainbow

Pastel

Fruit Punch

Tip: processing large data sets

2009-03-07T11:03:00.001-05:00

When deal with large data sets, ReadList is more efficient than Import.

From documentation:

If file is not already open for reading, ReadList opens it, then closes it when it is finished. If the file is already open, ReadList does not close it at the end.

ReadList[stream] reads from an open input stream, as returned by OpenRead.

So we can use OpenRead to open the file, then use ReadList to load large data sets piece by piece.

str = OpenRead[datafile]

(* the first number in our data file is the total number of points *)

ReadList[str, Number, 1]

{17581099}

Then we can read the real data in much smaller pieces every time, otherwise it may not have enough memory to run the processing algorithm.

Here is the example

(* read first 5 records *)

ReadList[str, Number, 7*5, RecordLists -> True] // MatrixForm

(* read the next 5 records *)

ReadList[str, Number, 7*5, RecordLists -> True] // MatrixForm

At the end, close the file by

Close[str]

Antialiasing in 3D graphics

2009-02-27T10:02:00.001-05:00

From documentation, this is a new feature in Mathematica 7.0:

“For 3D graphics, the operation of Antialiasing can depend on the particular graphics hardware you are using. Antialiasing is disabled unless Allow Antialiasing is set in the Preferences dialog.”

See the difference:

If you use a laptop with the integrated on-board video, antialiasing in 3D may not be supported at all.

A Bunny and many cuboids

2009-02-25T21:29:00.001-05:00

Counting coins with image processing

2009-02-22T20:12:00.001-05:00

Today I read a paper about automatically counting trees from satellite image, that makes me remember the coin-counting example from image processing class which I took many years ago.

I put some coins on the back of the mouse pad, snap a photo with the web camera.

img2 = Binarize[img]

img3 = ArrayPlot[MorphologicalComponents[img2]]

Now let’s do the counting.

cd = Tally[Flatten[MorphologicalComponents[img2]]]

{{0, 162643}, {1, 7600}, {2, 7633}, {3, 7369}, {4, 4309}, {5,
7679}, {6, 5868}, {7, 5883}, {8, 4082}}

(* pick the image region which is possible for the coins *)

cd2 = Select[cd[[All, 2]], # > 500 && # < 10000 &]

ListPlot[cd2, PlotStyle -> {PointSize[Large]}]

Clearly, there are three different groups of coins

Total[Map[Length, Sort[FindClusters[cd2]]] {5, 10, 25}]

130

Using $xxxx Mathematica to counting coins is kind of stupid itself, however, there are plenty of engineering applications actually do the very similar things in more complex situations.

Howto: Display 2D plot in 3D

2009-02-16T16:06:00.001-05:00

In GIS field, sometimes we like to stack several 2D plots together and display them inside a 3D box. For Mathematica, we need to define a function to convert a 2D plot into a 3D graphic object. I will use a small Geotiff as an example, you can download it here if you like to try the code.

(*get the ElevationRange, then import data *)

Import["smalldem.tif", {"Geotiff", "ElevationRange"}]

data = Import["smalldem.tif", {"Geotiff", "Data"}];

Then we a create the contour plot with 100 feet contours

c1=ListContourPlot[data,MaxPlotPoints->30,Contours->Function[Range[650,850,100]],ColorFunction->”DarkTerrain”,PlotRange->{640,850}]

Here the function we need to convert 2D plot into 3D

to3d[plot_,height_,opacity_]:=Module[{newplot}, newplot = First@Graphics[plot];newplot=N@newplot /. {x_?AtomQ,y_?AtomQ}->{x,y, height} ;
newplot /. GraphicsComplex[xx__]->{Opacity[opacity], GraphicsComplex[xx]}];

This function has three parameters: 2D plot, height, and opacity

Let’s create two more contour plots with 50 and 20 feet contours respectively. Then we can stack them together by setting them in different heights.

Show[{Graphics3D[to3d[c1,30,0.75]]}, Graphics3D[to3d[c2,20,0.75]], Graphics3D[to3d[c3,10,0.75]], Lighting->"Neutral", BoxRatios->{1,1,0.8},Axes->True]

Then we like to stack the original geotiff at the very bottom. This time we need to convert the raster into 3D. I use the example you can find in Listplot3D (check the section of “Neat Examples”).

r1=ReliefPlot[data,ColorFunction->colorf,ImagePadding->None, Frame->False, ImageSize->{800,800}];

pic = Reverse[ImageData[r1]];

bg=ListPlot3D[Table[1,{x,1,800,5},{y,1,800,5}],Mesh->None, VertexColors->{pic[[1;;800;;5,1;;800;;5]]},DataRange->{{1,800},{1,800}}, Lighting->"Neutral"]

The final product:

Here is the complete notebook.

Display large DEM with MaxPlotPoints

2009-02-13T16:40:00.001-05:00

MaxPlotPoints is a very useful option when works with large DEM.

dem=Import[tif,{"GeoTIFF","Data"}];
Dimensions[dem]
{2825, 5075}

Try ReliefPlot[dem], it takes a very long processing time in my not-up-to-date computer, on the other hand, with limited monitor resolutions, it doesn’t make sense to draw every point of this DEM.

Let’s first try:

ReliefPlot[dem, MaxPlotPoints –> 25]

Not quite good, the pixel is in rectangle shape rather than square, then with the small adjustment:

ReliefPlot[dem, MaxPlotPoints -> {50, 25}]

ReliefPlot[dem, MaxPlotPoints -> {1000, 500},
ColorFunction -> "GreenBrownTerrain"]

Click to see the large image, please ignore the left part, it is a quality issue from the DEM data source.

The same rule goes with ListPlot3D, if you want to get a 3D version.

ListPlot3D[dem, MaxPlotPoints -> {600, 300}, Mesh -> None,
Axes -> None, ColorFunction -> "GreenBrownTerrain",
BoxRatios -> {2, 1, 0.2}, ImageSize -> 1200]

Import water data from National Weather Service

2009-02-11T17:18:00.002-05:00

Besides real-time water data from USGS in the previous post, we like to import more data from National Weather Service, it has the forecast data on water stage.

The data is distributed in XML format.

NWSxml=Import[http://www.crh.noaa.gov/ahps2/xml/nori3_hgirg.xml];

(* function to extract data from xml *)

extractdata[xml_,group_]:=Module[{dataxml,groupdata},dataxml=Cases[NWSxml,XMLElement[group,_,_],Infinity];
groupdata=Transpose[Cases[dataxml,XMLElement["datum",{},{___,XMLElement[#,_,{x_}],___}]:>x,Infinity]&/@{"valid","primary","secondary"}];
groupdata /. {x_,y_,z_}:> {StringReplace[StringDrop[x,-9],{"T"->" "}],ToExpression[y],ToExpression[z]*1000}];

(* extract data*)

observeddata = extractdata[NWSxml, "observed"];
forecastdata = extractdata[NWSxml, "forecast"];

(* produce the graph *)

Flood stage is at 11.0 feet, it seems the white river may be in the flood stage in next several days.

Import USGS national real-time water data

2009-02-10T22:54:00.003-05:00

The built-in real-time weather data in Mathematica 7.0 is wonderful. You may be interested in importing more related data.

USGS has web-interface for Daily Streamflow Conditions

Here is the simple steps from importing Gauge height and Discharge data

USGS 03353800 White Lick Creek at Mooresville, IN

(*import html as plain text *)

rawdata=Import[“http://waterd ata.usgs.gov/nwis/uv?cb_00060=on&cb_00065=on&format=html&period=7&site_no=03353800”];

(* extract table information *)

rawdata=StringReplace[rawdata,","->""];
result=StringCases[rawdata,x:DatePattern[{"Month","/","Day","/","Year"," ","Hour",":","Minute"}]~~Whitespace~~y:NumberString~~Whitespace~~z:NumberString:>{x,ToExpression[y],ToExpression[z]}];

(*draw graph*)

DateListPlot[result[[All,1;;2]],PlotStyle->{Blue, PointSize[Medium]},PlotLabel->"Site: 03353800",Joined->True,Filling->Bottom,FrameLabel->{"Time","Gage Height (feet)"},DateFunction:>(DateList[{#,{"Month","/","Day","/","Year"," ","Hour",":","Minute"}}]&),DateTicksFormat->{"MonthShort","/","Day"," ","Hour",":","Minute"}]

By the way, the code is tested in Mathematica 6.0.

Create the map with a image background

2009-01-05T00:14:00.001-05:00

Here is something you can do with Mathematica 7.0 Image Processing

Here is a good picture of American Flag

Then create a mask for it and add them together:

mask= Graphics[{Black, CountryData["UnitedStates", {"Polygon", "Mercator"}]}, ImageSize -> {500, 335}]

flagImage = ImageAdd[img, mask]

Then create another new image with the boundary:

bgk = Graphics[{FaceForm[], EdgeForm[{Black}], CountryData["UnitedStates", {"Polygon", "Mercator"}]}, ImageSize -> {500, 335}]

ImageMultiply[bgk, flagImage]

How to grab photos from Flickr

2008-12-31T22:07:00.001-05:00

Happy new year 2009!

I have a post about make a Flickr image wall with Mathematica

First, you need get some basic knowledge about Flickr API, Photo Source Url

The coding is very simple, just play with it, then you can modify it to try something new.

catxml = Import[
"http://www.flickr.com/services/rest/?method=flickr.groups.pools.getphotos&api_key=your_API_key&group_id=75601600@N00&per_page=500&page=1&safe_search=1", "XML"];

This line will do the query from Flickr group Somebody else's cat, return 500 images’ meta data

cats = Cases[catxml, XMLElement["photo", id_, {}] :> id, Infinity];
ids = #[[All, 2]] & /@ cats;
ids = RandomSample[ids, 100];

This three lines will extract 100 random photo ids

cs="http://farm"<>#[[5]]<>".static.flickr.com/"<>#[[4]]<>"/"<>#[[1]]<>"_"<>#[[3]]<>"_s.jpg" &/@ids;
imgs=Import[#] & /@ cs;

Then get the real url to each photo and import them. Here we use small size (75 by 75) “_s.jpg”

ImageAssemble[Partition[imgs, 10]]

Then put them together into one image

So we finish it just in several lines.

How to make a Ternary Plot

2008-12-11T13:51:00.004-05:00

For the reader who are curious about the post on Ternary Plot

The code is very simple, here is the mathematica notebook, and the pdf version if the notebook somehow doesn't work.

Sorry for the late reply, I seldom check the comments for the older posts.

Tip: create 3D Map

2008-12-07T13:02:00.001-05:00

Here is a simple explanation on how to create 3D map in the previous post.

The key is to use Extrude function (I copied it somewhere):

extrude[pts_, h_] := Module[{vb, vt, len = Length[pts], nh},
If[! NumericQ[h], nh = 0., nh = N@h];
vb = Table[{pts[[i, 1]], pts[[i, 2]], 0}, {i, len}];
vt = Table[{pts[[i, 1]], pts[[i, 2]], nh}, {i, len}];
GraphicsComplex[
   Join[vb, vt], {Polygon[Range[len]],
    Polygon[Append[
      Table[{i, i + 1, len + i + 1, len + i}, {i, len - 1}], {len, 1,
       len + 1, 2 len}]], Polygon[Range[len + 1, 2 len]]}]]

Let's create a map to show per-capita oil consumption in South America.

(* get the data *)

(* Falkland Islands is removed *)

mapdata = {First@CountryData[#, "Polygon"],
     CountryData[#, "OilConsumption"]/
      CountryData[#, "Population"]} & /@
   DeleteCases[CountryData["SouthAmerica"], "FalklandIslands"];
mapdata[[All, 2]] = Rescale[mapdata[[All, 2]]];

(* 3D map *)

Graphics3D[{EdgeForm[Darker[ColorData["Rainbow"][#[[2]]]]],
FaceForm[ColorData["Rainbow"][#[[2]]]],
extrude[#[[1, 1]], #[[2]]]} & /@ mapdata, Boxed -> False,
Lighting -> "Neutral", BoxRatios -> {1, 1, 0.2}, ImageSize -> 800]

The visual quality of 3D map is low, however, it is good enough for the classroom.

Mathematica 7: import ESRI Shape file

2008-12-03T19:57:00.001-05:00

Finally, ESRI Shape file is supported in Mathematica 7: Import["dir/file.shp"]

However, I feel it is necessary to explain it in a little more detail for users who are not that familiar with Mathematica.

Let's import a shape file from Indiana GIS Atlas which contains county boundaries.

(*import data*)
data=First@Import["http://129.79.145.7/arcims/statewide_mxd/downloads/LANDSURVEY_COUNTY_POLY_IN.zip", Data"];
(*check what is inside the data*)
data[[All, 1]]
{"LayerName", "Geometry", "Labels", "LabeledData"}
(* check the field name in attribute table *)
"Labels" /. data
{"AREA", "PERIMETER", "NAME_U", "NAME_L", "NCAPC"}
(* attribute data is stored in "LabeldData" *)
area = "AREA" /. ("LabeledData" /. data);
area = Rescale[area];
(* get all the geometry *)
geometry = "Geometry" /. data;
(* put two data together *)
mapdata = Transpose[{geometry, area}];
(* generate the thematic map *)
Graphics[{EdgeForm[{Thick, Gray}], FaceForm[ColorData["Rainbow"][#[[2]]]], #[[1]]} & /@ mapdata, Frame -> False, ImageSize -> 400]

With some extra efforts, a 3D version is generated.

Extract font outlines

2008-11-24T23:34:00.001-05:00

I find this tip from the talk: Doing Physics with Mathematica 6 by Michael Trott

GCText[text_,tsize_]:= ImportString[ExportString[Cell[TextData[StyleBox[
text, FontWeight->"Bold",FontFamily->"Helvetica"]], "Text", tsize], "PDF"]][[1]]

outlines = (First /@ Cases[GCText["ABC", 20], _Polygon, Infinity]);

Graphics3D[{FaceForm[],EdgeForm[Black],Polygon[Map[Append[#, 10]&, outlines, {-2}]]}];

for fun:

QHull for Mathematica

2008-10-05T22:54:00.001-04:00

mPower.m is a Mathematica package that interfaces with qhull binaries.

Installation guide is not for windows platform.

Here is the quick way to get it running on Windows Vista with Mathematica 6.0.

1. download mPower.m, and put it into the Applications subdirectory of $UserBaseDirectory

In[1]:= $UserBaseDirectory

Out[1]= C:\Users\...\AppData\Roaming\Mathematica

2. download qhull for windows, put it into the folder C:\qhull

3. open mPower.m, modify the following two lines,

$QHULL=ToFileName[{$UserBaseDirectory, "Applications", "qhull", "src"}]

to: $QHULL="C:\\qhull\\"

qhullFiles={"qconvex", "qdelaunay", "qhalf","qhull", "qvoronoi"};

to: qhullFiles={"qconvex.exe", "qdelaunay.exe", "qhalf.exe","qhull.exe", \
"qvoronoi.exe"};

4. then testing, ignore the warning on regtet and pwrvtx, we only need qhull

<< "mPower.m"

points = RandomReal[1, {10, 3}];

ch = convexHull[points, convexHullFormat -> {facetNormals -> True, facets -> True}]

Thanks to the mPower developers

Financial Data: AAPL vs MSFT

2008-08-07T13:37:00.002-04:00

Technologically speaking, Apple and Microsoft are two quite different companies in many ways. We all know that in recently several years, Apple stock is very hot. Here is the stock price: APPL vs MSFT since 2004, the legend is not necessary for this graph.

DateListPlot[{FinancialData["AAPL",{2004}], FinancialData["MSFT",{2004}]}, Joined->True, Filling->Bottom]

Google Insights for Search gives another way to present it.

Google search volume index for "AAPL" since 2004

Google search volume index for "MSFT" since 2004

How to extract points from a plot

2008-08-06T14:49:00.001-04:00

pts=Cases[theplot, x_GraphicsType:>First@x, Infinity];

Example 1:

g=Plot[Sin[x],{x,0,2Pi}];
pts=Cases[g,x_Line:>First@x,Infinity];
ListPlot[pts]

Example 2:

g=Plot3D[Cos[x]+Cos[y],{x,0,4Pi},{y,0,4Pi}];
pts=Cases[g,x_GraphicsComplex:>First@x,Infinity];
ListPointPlot3D[pts]

Wolfram Screencast Gallery

2008-07-16T09:56:00.001-04:00

There are some very impressive presentations in Screencast Gallery.

I watched Spinthariscope by Theodore Gray. This simple shape editor is particularly interesting.

Fun with Gamepad

2008-07-11T09:41:00.003-04:00

We have a cheap Logitech gamepad (less than 10$ from Target).

ControllerInformation[]

Controller Device 3: Logitech(R) Precision(TM) Gamepad
Raw Product Name "Logitech(R) Precision(TM) Gamepad"
Raw Product ID 49690
Device Type Windows DirectInput Device
Raw Controller Type Joystick
Mathematica Controls 22 controls
Raw Controls 12 controls

Using gamepad to rotate 3D graphics is very smooth, and can jump to top view from x, y, z directions when you hold down different buttons. It seems not supported to zoom a 3D graphic with gamepad.

Gamepad & Device Interface for the details.

Google Earth is another application to play nicely with gamepad.

Update: See here for enabling zoom function with gamepad. Thanks.

Standard Extra Packages

2008-07-09T11:22:00.001-04:00

See here Mathematica 6.0 Standard Extra Packages

It gives out a very helpful overview on packages not loaded by default.

T-Thirt design: Pi with rainbow color scheme

2008-06-15T17:10:00.003-04:00

Functions used: Rasterize, Grid

Something hidden

2008-05-08T10:57:00.004-04:00

No idea? Look at the one:

Rasterize

2008-05-04T18:13:00.001-04:00

Self random walk on sphere

2008-05-03T22:42:00.007-04:00

Self random on sphere lattice (longitude, latitude).
(red -- starting point, green -- ending point)

Draw on a sphere:Self avoid version:

The longest horizontal lines indicate the jumps between -180 and 180.

Opacity

2008-04-27T15:35:00.002-04:00

Nationality Name

2008-04-21T11:49:00.005-04:00

For most counties, the nationality name usually in form of country name + suffix: -an, -ian, and -ese.
China->Chinese, Japan->Japanese, are all the east Asian the "ese" people? Let's find out with Mathematica: CountryData["Country", "NationalityName"].

From the above map, you can see it is not true at all, e.g. Korea->Korean, Thailand->Thai. Portugal->Portuguese is only "ese" country in Europe.

One thing I hear quite often is that -ese means smaller and less power as opposite to -an and -ian, -ese is used by western colonist in an insulted way. I don't know how much historic truth in this claim. Google search gives out the following:
-ian -an from Latin –ianus, meaning "native of", "relating to", or "belonging to"
-ese from the Latin -ensis, meaning "originating in"
However, it is hard to know why some countries are with "-ese", others are not.

Nonsense: Something here doesn't belong

2008-04-20T15:06:00.004-04:00

Nonsense

2008-04-10T22:02:00.003-04:00

Rejected logo design in the summer, 1999.

Colorful Words

2008-04-10T21:22:00.003-04:00

Rasterize, then coloring

Minimum Bounding Rectangle of a Point Set

2008-04-07T23:16:00.006-04:00

Give a 2-D point set, we need to find the minimum bounding rectangle in term of area.

I read it somewhere that the minimum rectangle must always contain at least one edge of the convex hull, however, I can't find a citation now. So the algorithm can be constructed in the following way: first construct the convex hull, then rotate each edges for the convex hull to the position parallel to x-axes, then calculate area of bounding rectangle; find the minimum one of these rotated rectangles and rotate it back to normal.

Nonsense

2008-04-04T13:42:00.003-04:00

The early sign of scientist self-destruction nature.

mushroom is from http://www.wpclipart.com/plants/mushroom/mushroom.png

Google Chart API

2008-04-04T10:25:00.003-04:00

The Google Chart API lets you dynamically generate charts.

The basic format:
http://chart.apis.google.com/chart?<parameters>

It is the most straightforward solution for web application to display dynamic graphs. I am wondering if Wolfram can come up something similar, not the webMathematica, just a simple free image generator will be great, it is definitely a good way to get more publicity.

Tips: Customizing Graphplot

2008-04-03T10:22:00.004-04:00

For a graph g, we like draw leaves in different style from nodes.

(* find leaves *)
leaves = Complement[g[[All, 2]], g[[All, 1]]];
TreePlot[g, VertexLabeling -> False, PlotStyle -> {Black}, VertexRenderingFunction -> (If[MemberQ[leaves, #2], {FaceForm[LightGray], EdgeForm[Black], Disk[#1, 0.15], Text[#2, #1]}, {FaceForm[White], EdgeForm[Black], Disk[#1, 0.15]}] &), AspectRatio -> 0.3]

The tip is in VertexRenderingFunction, Text[#2, #1], #2 means the label, #1 is coordinates.

Learning regular expression

2008-03-22T16:01:00.003-04:00

I am learning Regular Expression right now. Here is the example to pull the information on monthly posting activities of Mathematica Discussion Group.

url = "http://groups.google.com/group/comp.soft-sys.math.mathematica/about";
reg = "(?<=month/)\\d+-\\d+\">\\d+";
data = Import[url, "Source"];
ex = StringCases[data, RegularExpression[reg]];
ds = ToExpression[StringSplit[ex, {";", "-", "\">"}]];

Then we can use DataListPlot to produce a nice graph.

Understanding Finanical Market: Arbitrage

2008-03-12T12:26:00.006-04:00

From Wikipedia, Arbitrage is the practice of taking advantage of a price differential between two or more markets, the profit being the difference between the markets prices.

Fro example, China Petro & Chem is traded in both US and HK stock markets. The price between these two markets are different, it seems US price follows the movements of HK's in recent several months.

After converting HK share price to US dollars, the price differences can be shown:
DateListPlot[{ussnp, hksnp}, Joined -> True, Filling -> {1 -> {2}}, FillingStyle -> {Green, Red}, DateTicksFormat -> {"Month", "/", "YearShort"}]
The trick is to use Filling -> {1 -> {2}}, FillingStyle -> {Green, Red}, Green indicates HK price is higher than US price, the red means the opposite.

Tip: Orthographic view of 3D graphics

2008-03-03T14:17:00.010-05:00

The orthographic view of the 3D plot can be generated by using ViewPoint option. It comes very handy in some situations.

data = RandomReal[1, {20, 3}];
ListPlot3D[data, InterpolationOrder -> 0, ColorFunction -> "Rainbow", Mesh -> None, Boxed -> True, Axes -> None]
Add option ViewPoint->{0,0,Infinity} and draw it again.

Then we can put together a nice bounded Voronoi diagram. It saves the trouble to use BoundedDiagram from Computational Geometry Package.

Min Max problem

2008-03-02T09:21:00.003-05:00

There is a discussion in Matlab blog: Should min and max marry?
The question is what the overhead is in calling min and max separately, instead of scanning through the data once, looking for both the lowest and highest values, and is it worth having a combined function for min and max from a speed point of view? It seems worthy it, if you work with large data frequently.

I am not sure about how Mathematica handle Min/Max internally. The one thing I can confirm is that Min/Max operations on large random 1 dimensional array are several times faster than the ones in Matlab in my computer.

Tip on Mathematica Programming

2008-02-25T10:37:00.002-05:00

There are two interesting presentations from 2007 Wolfram Technology Conference.

Efficient Mathematica Programming: tips from MathGroup, very helpful.
A New Mathematica Programming Style: prefix and postfix-pure function is disscussed as an alternative to the traditional matchfix-dominated style

Google Static Maps API

2008-02-25T09:56:00.007-05:00

Google finally release Google Static Maps API, which allows you to generate the maps using a regular URL (ala REST) along with parameters specifying location, size, etc and it returns a unique GIF image with that map.
However, there are several drawbacks, for example, 512x512 is the largest image size allowed and geocoding is not integrated. Even though, it is better than nothing.

The following shows the API example:
Top 10 "Mathematica" cites from Google Trends.

For fun, top 10 Mathematica cities vs top 10 Matlab cities
Red "M"--> Mathematica City, Blue "T" --> Matlab City, Green "B" --> City with both Mathematica and Matlab.

Update: static maps seems use the old data
Noticed the problem with "Trade Center Dr",
Check this place on Google Map.

Play with Images: Letter Image

2008-02-17T21:37:00.005-05:00

Spell with Flickr is interesting. I modified the Flickr image wall code.
Here is the result:

Import High Resolution DEM (2)

2008-02-08T09:23:00.000-05:00

We are talking about more general case here: DEM is processed in ArcGIS (e.g. Mosaic, Extract by Mask, etc.). And it is ready for computation in Mathematica. The DEM can be exported as Float32 format, imported with the procedure in the previous post.
In this post, the DEM is exported as ARCGrid ASCII format, the whole DEM is in one text file.

data = Import["DEM.txt", {"Lines", All}];
(* the first lines are metadata *)
meta = data[[1 ;; 6]]

{"ncols 1703", "nrows 1200", "xllcorner
-86.528233463385", "yllcorner 41.430023494444", "cellsize
0.00027777777779994", "NODATA_value -9999"}

(* drop the metadata *)
data = Drop[data, 6];
(* convert string to data *)
data = Map[ToExpression, StringSplit[#]] & /@ data;

The next step is that we will make the output of ReliefPlot match the other geographic data.

1. Deal with NOData: PlotRange->{zMin,zMax}

2. Put the image into right geographic zone: DataRange-> {{Lon_Min, Lon_Max}, {Lat_Min, Lat_Max}}

ReliefPlot[data, DataReversed->True, DataRange-> {{-86.5282, -86.0552}, {41.4300, 41.7634}}, PlotRange -> {199.0, 283.0}]

Then tested it with the map created in another post.

Note: ReliefPlot needs an option: NoData->{}

Import High Resolution DEM (1)

2008-02-08T08:45:00.001-05:00

The best place to get high resolution DEM is from USGS National Map Seamless Server.

It publishes National Elevation Dataset (NED) in 1 arc second(~30m resolution), 1/3 arc second(~10m resolution) and 1/9 arc second(~3m resolution).

The default file format is ESRI ArcGRID. Several other non-proprietary formats are available. If you don't want to further process the DEM in ArcGIS, GridFloat format is the right choice. GridFloat is "a simple binary raster format (floating point data). There is an accompanying ASCII header file that provides file size information (number of rows and columns). The data are stored in row major order (all the data for row 1, followed by all the data for row 2, etc.)". See National Elevation Data (NED) FAQ.

Follow this tutorial for downloading procedure.

After downloading, unzip the file, .flt is the data for DEM, .hdr is the header ASCII file.

Open header file (.hdr)

ncols 650
nrows 395
xllcorner -86.428333330255
yllcorner 41.639166666159
cellsize 0.00027777777779995
NODATA_value -9999
byteorder LSBFIRST

We can see this DEM has 395 row and 650 col.

In Mathematica:

data30 = Import["ned30.flt", {"Binary", "Real32"}];
(* each row has 650 col. *)
data30 = Partition[data30, 650];
ReliefPlot[data30, DataReversed -> True, ColorFunction -> "Rainbow"]

The following examples show the same area in 1 arc second (30m) and 1/3 arc second (10m) DEMs.

30m

10m

Select Point Data inside Geographical Boundary

2008-02-07T11:08:00.000-05:00

Here is an example on how to select point data inside the geographical boundary by using point inside polygon function.

(* all the cities {name, coordinates} in Indiana *)
allcities = {First[#], CityData[#, "Coordinates"]} & /@ CityData[{All, "Indiana", "UnitedStates"}];

We like to select the cities inside St. Joseph county.
Boundary of St. Joseph county is imported as Polygon. See details on importing GIS data.

In Mahtematica-users Wiki, there is an discussion on Point Inside Polygon. We choose the first solution, since it works on non-convex polygon.

(* select cities by Point inside Polygon function *)
Select[allcities, PLSPolygon[countyboundary, #[[2]]] &]

Fun with Financial Data: Currency Exchange Rates

2008-02-05T11:03:00.000-05:00

Someone posted some interesting photos from his trip to Japan in this January. One is delicious sushi sold at the roadside in somewhere, Kyoto, Japan.

In the top row, the one marked for 900 Yen looks exactly like the ones sold here in greocery store. I am wondering the price in US dollars.

900/FinancialData["USD/JPY"] gives $8.42, it is around 60% more expensive than the ones you can get here, however, it must taste much better.

One interesting thing is how much we can get if we convert 1 USD dollar to Japan Yen through the third currency.

FinancialData[{"USD",#}] / FinancialData[{"JPY",#}] &/@ {"HKD","EUR","GBP","CAD"}
{106.842, 106.625, 106., 106.947}

Direct exchange rate is 106.825.

USD->GBP->JPY: lost near 1% USD->CAD->JPY: gain tiny 0.1%

Is there a simple way to find a currency exchange chain "USD->..->... ->JPY" to maximize the gain over "USD->JPY"?

Fun with MatrixPlot: Pixel Art

2008-02-04T12:55:00.000-05:00

MatrixPlot is great for pixel art.

jp=MatrixPlot[pixelart, Frame->None, ColorRules->{0->White}, ColorFunction->#, ImageSize->90, Mesh->All, MeshStyle->LightGray] &/@ Append[ColorData["Gradients"], "Monochrome"]; GraphicsGrid[Partition[RandomSample[jp],8],Spacings->0]

ColorRules can be used for the better result, e.g. ColorRules -> {0 -> White, 1 -> Yellow, 2 -> Orange, 3 -> Pink, 4 -> Red}

Tip: Import 3D Model

2008-01-31T21:25:00.000-05:00

For whatever reason, some 3D models have to be imported in Mathematica.
The 3D model is stored as OBJ file. When it is imported, it loses the texture information.

(* set the color *)
roof = Import["greenroof.obj"]
roof = roof /. {l_Line :> {Thick,Black, l}, p_Polygon :> {FaceForm[Darker[Green]], p}}

Here the whole model:

Tip: Mathematica version of Maltab "Find" Function

2008-01-29T09:56:00.000-05:00

I must say that if you work with Matlab everyday, it is really frustrated with Mathematica in some points.

Find is one of the most frequently used functions in Matlab.
Find(array>5.0): returns indices of the elements that are greater than 5.0.

In Mathematica, the corresponding function is Position

Position[array, #>5.0&] or Position[array, _>5.0] just doesn't work.

We "simple-minded" Matlab user have seldom thought about putting constraints on patterns

Then, this will do the job

Position[array, x : _ /; x > 5.0]

Position[array, x:_ /; testfunc[x]] will do better.

Tip: ListPlot3D and InterpolationOrder

2008-01-28T12:05:00.000-05:00

Not exactly a tip.

10,000 random points
Regular grid: data = Table[RandomReal[], {100}, {100}];
Irregular grid: data = RandomReal[1, {10000, 3}];

Table[ListPlot3D[data, InterpolationOrder -> n, Mesh -> None]; // Timing, {n, {0, 1, 2, 3, 4}}]

Result (in seconds):
Interplotion Order: 0, 1, 2, 3, 4
Regular Gird: 0.771, 3.415, 3.435, 3.525, 3.565
Irregular Grid:6.349, 4.871, 4.767, 4.777, 4.737

Difference in zero-order interpolation is huge.

Reference: InterpolationOrder

Tip: Counting Files by Date

2008-01-24T15:33:00.000-05:00

We like to know the number of the certain special type file processed month by month.
The following code will do the job:

(* all the files under one base direcotry *)
SetDirectory[basedirectory]
(* find all the files in all the subdirectories*)
files = FileNames["*.byt", {"*"}, Infinity];
(* return {year, month} for each file *)
data = Take[FileDate[basedirectory <> "\\" <> #], 2] & /@ files;
(* sort then count *)
cs = Tally[Sort[data]];
(* draw the graph *)
DateListPlot[cs, Filling -> Bottom, FillingStyle -> Blue, PlotStyle -> {Red, PointSize[Large]}, PlotRange -> {{{2006, 1}, {2008, 2}}, Automatic}, DateTicksFormat -> {"MonthShort", "/", "YearShort"}]

Fun with Map: Simple Geocoding with Yahoo Map API

2008-01-22T11:53:00.001-05:00

Yahoo Map API has a simple REST interface which works well with Mathematica.
Here is an example to show a simple Geocoding function.

Problem: Give an address, return (longitude, latitude)

The basic query format:
http://local.yahooapis.com/MapsService/V1/geocode?appid=YahooDemo
&street=?&city=?&state=?

Let's try it with the address of Wolfram Research Inc.

geocoding[{"100 Trade Center Drive", "Champaign", "IL"}]
{-88.244868, 40.097855}

Is the result right? Let's check it out with Map Image API.

image[geocoding[{"100 Trade Center Drive", "Champaign", "IL"}]]

Notice:
1. I use "YahooDemo" as AppID, if you want to use this function frequently, please apply your own ID, it is free.
2. The Geocoding API is limited to 5, 000 queries per IP per day, Google Map has the similar limitation, so if you need geocoding a large amount of addresses, please use offline geocoding application.

reference: Yahoo Map API, Map Image API, Geocoding API

update:
Google version: see 1th comment , Thanks very much.

Tip: Programming Style

2008-01-20T22:42:00.000-05:00

Read it from the company newsletter sometime ago:

Given a list of pairs of numbers, return a list consisting of the sum of each pair.
pairs = {{58, 96}, {85, 22}, {100, 69}, {5, 37}, {32, 64}, {41, 86}, {14, 0}, {79, 22}, {55, 36}, {86, 39}, {11, 15}};

(* Style I*)
result = Table[Null, {Length[pairs]}];
Do[result[[k]] = pairs[[k, 1]] + pairs[[k, 2]], {k, 1, Length[pairs]}]
result
(* Style II *)
Table[pairs[[k, 1]] + pairs[[k, 2]], {k, Length[pairs]}]
(* Style III *)
Apply[Plus, pairs, {1}]
(* Style IV *)
Map[Total, pairs]
(* Style V *)
pairs /. {p_, q_} -> p + q
(* Style VI *)
pairs[[All, 1]] + pairs[[All, 2]]

I work with Maltab most of the time, it is useful to read this kind of document occasionally.

Details on Importing GIS Data

2008-01-18T10:08:00.000-05:00

If you live US, there are plenty of GIS data published by different government agencies, it is quite possible that the county where you live has its own GIS department now. Most of the vector data come in ESRI shape file format. There are several different ways to get shape file into Mathematica. One simple way is explained here, it uses tools based on shapelib (open source tools to read/write shape file).

Step 1. Find the data you want
Just go to the agency may have the data you are interested. In this example, I will use the Indiana county boundary data I have.
There are some free/open source GIS software you can use to check the data: Quantum GIS, ArcGIS Explorer. (No offense here, I just assume you don't much about GIS.)

Step 2. Output shaple file to data
We use shp2text, and you can find more information on its website.
In command line mode:
shp2text --gpx indiana.shp
the output is the structure of the associated data in DBF file. We like to have county name, it is in Field 1.
Field 0: Type=Integer, Title=`FIPS', Width=5, Decimals=0
Field 1: Type=String, Title=`COUNTY', Width=16, Decimals=0
Field 2: Type=String, Title=`DOQ_DISK', Width=16, Decimals=0
Field 3: Type=String, Title=`FIPS_STRG', Width=16, Decimals=0
Field 4: Type=String, Title=`ATLAS_CODE', Width=16, Decimals=0
Field 5: Type=Double, Title=`AREA_FT2', Width=19, Decimals=3
Field 6: Type=Double, Title=`AREA_ACRES', Width=19, Decimals=3
Field 7: Type=Double, Title=`PERIMETER_', Width=19, Decimals=3

shp2text --gpx indiana.shp 1 0 > output.xml

There will dump the shape file into xml file, of course, if you can dump it into plain text without "--gpx", and reformat the plain text to cvs format. I personally prefer to XML format.

see the whole document here.

Step 3. Process XML file in Mathematica
I listed the code step by step here in case you are not familiar with XML.

dataxml = Import["output.xml"];
(* extract county names *)
countyname = Drop[Cases[dataxml, XMLElement["name", __, {w_}] :> w, Infinity], 1];
(* extract {lat, lon} for each county boundary *)
wholes=Cases[dataxml, XMLElement["rte", __], Infinity] ;
counties = (Cases[#, XMLElement["rtept", __], Infinity] & ) /@ wholes;
dataLatLon = counties /. (XMLElement["rtept", latLonRules_, {}] :> { {"lat", "lon"} /. latLonRules });
(* tidy up a little, of course, you can combine them into one line code*)
test = Flatten[#, 1] & /@ dataLatLon;
test = Map[ToExpression, test];
test = Map[Reverse, test, 2];
(* then draw it *)
Graphics[{EdgeForm[Gray], FaceForm[LightGreen], Polygon[test]}]

Here is the example that uses the data imported as the background map.
Point data is from CityData, and it shows all cities with population > 1000.

Importing point, line type of data is in the similar way.

Tip: Efficient Plot Large Point Data Set

2008-01-18T09:33:00.000-05:00

It is from the discussion from math group.
Point, Line and Polygon now accept multi-coordinate list that can be processed and rendered more quickly by the Mathematica front end than the equivalent individual primitives.

data = Table[{Random[], Random[]}, {100000}];
Graphics[{PointSize[0.002], Point[data]}, AspectRatio -> Automatic]
Graphics[{PointSize[0.002], Point[#] & /@ data}, AspectRatio -> Automatic]

More information: Efficient Representation of Many Primitives

Old Time Fun: Ternary Plot

2008-01-17T14:40:00.001-05:00

Here is the ternary plot of the 140 named colors supported by modern browsers. Actually, only 139 colors are drawn in the graph. There is no place for Black (RGBColor[0,0,0]), the center is White, however, you can't see it on the white background.

Old Time Fun: Binary Image

2008-01-15T14:25:00.001-05:00

Can you find out the relationship between these two images?

Old Time Fun: Hilbert Curve

2008-01-14T11:15:00.000-05:00

At the time I started learning Mathematica (maybe version 3.0), it is almost impossible for any student not to plot space filling curves or other fancy fractals. Except for eye candy, it is quite useless for most of us. Today, it can't even be counted as eye candy any more.
Here is mine, Hilbert Space Filling Curve, which I found in pile of old files, the only thing that warms my heart a bit is that I come up with 4 dimensional version. You probably don't see 4D ones quite often. Of course, it has to be projected into 3D in some way for displaying.

Play with Images: Sky Color

2008-01-11T15:04:00.001-05:00

I grabbed some images tagged with "sky" from Flickr with yesterday's image wall code. One thing interests me that if there is a simple way to tell the color of sky. I am not a digital image person. Let's us see what a layman can do abut it during the lunch break.
The most simple idea is that averaging the color of the upper 3/4 part of the image. After tried with several images, the result of averaged RGB value are mostly some gray colors. The HSB color system seems to be a much better choice.

The averaged H(hue) shall do the job. Here is the some examples with Quartiles[h](1/4, 1/2, 3/4) and Mean[h].

Median probably works better than Mean does for identifying the main sky color, and the color changes of different quartiles may be used as the indicator about the cloudy sky. I don't have time to go further.
note: RGBtoHSV is part of the image processing package.

Play with images: Flickr image wall

2008-01-10T11:28:00.000-05:00

I come across one old python script to download the images from flickr. It will be cool if it can be done in Mathematica. The goal at this lunchtime is to grab some thumbnails from group Field Guide: Birds of the World to make a small image wall.

Here is the output

Formated XML output

Image wall

Images are placed in random orders.

Achieved in 9 lines of code, not too bad.

reference: Flickr API, Photo Source Url

Created a larger version (100 photos).

Fun with Flags: who have star in their flags?

2008-01-08T10:00:00.000-05:00

How many countries have at least one star in their flags?

(* get flag descriptions of all the countries *)
desc = {#, CountryData[#, "FlagDescription"] /. _Missing -> ""} & /@ CountryData[];
(* select countries which have the "star" *)
pc = Select[desc, StringCount[#[[2]], "star"] > 0 &];
(* take the country name out *)
cc = pc[[All,1]];
(* take all the flags *)
flags = CountryData[#, "Flag"] & /@ cc;
(* show all the flags *)
GraphicsGrid[Partition[flags, {10}], ImageSize -> 800, Frame -> All]

There are total 80 countries and areas that have at least one star in their flags.

Play with Numbers: January Financial Market Performance

2008-01-04T12:58:00.000-05:00

I am getting curious about the stock market during the lunch time. I like to check the January performance in the past.

(* define the function to get the January data *)
cc[stock_, year_] := FinancialData[stock, {{year, 1, 1}, {year, 1, 31}}]
(* plot multiyear graphs*)
Table[DateListPlot[cc["^DJI", year], PlotStyle -> {If[First[Last[#][[2]] - First[#] [[2]] >= 0 & /@ {cc["^DJI", year]}], Green, Red]}, Joined -> True, PlotLabel -> year], {year, 1991, 2006}]

Here is DOW and NASDAQ January performance from 1991 to 2006. By comparing the first and last trading day, if it gains, the line color is green, otherwise, it is red.

DOW January Performance (1991~2006)

NASDAQ January Performance (1991~2006)

You see, Mathematica does come in handy sometime!

Overlay GraphPlot on maps

2008-01-03T10:34:00.001-05:00

GraphPlot is one of the most interesting and powerful function in Mathematica. Here comes up a simple demonstration to show how to overlay the Graph on top of the map. The key is to use VertexCoordinateRules option to place the each vertex in the right place.

The following graph is the adjacency graph for countries in South America.
Using the location of capital cities as vertex coordinates seems a good choice. So combining CountryData and CityData, the location list of {country->{lon, lat of capital city}} is created.

{{"Argentina" -> {-58.37, -34.61}, "Bolivia" -> {-68.15, -16.5}, "Brazil" -> {-47.91, -15.78}, "Chile" -> {-70.64, -33.46}, "Colombia" -> {-74.09, 4.63}, "Ecuador" -> {-78.5, -0.19}, "FrenchGuiana" -> {-52.34, 4.92}, "Guyana" -> {-58.16, 6.79}, "Paraguay" -> {-57.63, -25.3}, "Peru" -> {-77.05, -12.07}, "Suriname" -> {-55.2, 5.85}, "Uruguay" -> {-56.17, -34.87}, "Venezuela" -> {-66.93, 10.54}}

Redraw the graph with VertexCoordinateRules -> locationlist

Then show it with the map

It is kind of silly in certain way to overlap the adjacency graph on top of the map. However, the graph delivers the information more clearly, you can spot it immediately that Brazil has the largest number of neighbors in South America.

Check weather with Mathematica

2008-01-01T17:19:00.000-05:00

wxml = Import["http://www.weather.gov/data/current_obs/KBMG.xml", "XML"];

First[Cases[wxml, XMLElement[#, _, {w_}] :> w, Infinity]] & /@ {"weather", "observation_time"}

{"Light Snow", "Last Updated on Jan 1, 4:53 pm EST"}

Happy New Year 2008: in security image way

2007-12-31T20:50:00.001-05:00

ArrayPlot[ReplacePart[Rasterize[Framed["HAPPY NEW YEAR 2008"], "RGBColor", ImageSize->300, Background->LightYellow], Table[RandomInteger[{1,350},2] -> Hue[RandomReal[1]], {30000}]], Frame->False]