Dec 21, 2010

Happy Holiday!

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

Nov 15, 2010

Color-coded contour lines with Mathematica

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.

Nov 8, 2010

Customizing DateListPlot with PlotMarkers

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.

Nov 1, 2010

Plotting on Google static map

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:


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

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.

Oct 26, 2010

Processing comic images for e-reader

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:

    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.

Sep 29, 2010

Tips: Import MATLAB MAT-files

You want to import the “.mat” file:


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)

Aug 25, 2010

Work with Web Map Service (WMS)

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, 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.,-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

Aug 20, 2010

Visualize weather pattern with random walk

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.

Apr 19, 2010

Play with Bibliographic Data

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:

  <style face="normal" font="default" size="100%">Paul, N.</style>
  <style face="normal" font="default" size="100%">Cao, B.</style>
… </authors>

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.

Mar 25, 2010

Extract elevation data with Google Elevation Service

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

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:

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

ToExpression@StringCases[json, NumberString]


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.


Jan 26, 2010

More on Heatmap

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



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.


Jan 12, 2010

Charting time series as calendar heat maps

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.