Mar 26, 2014

Geospatial function: Point in Ploygon Test

There is an excellent post on point in polygon test from Mathematica Stackexchange.

Point in polygon test is one of the most useful functions when processing Geospatial data. Unfortunately, there is no official built-in function from Mathematica yet. Here is an example of showing usefulness of this function.
We have some data on Greenland icesheet thickness in format of {lon,lat, thickness}.


Our goal is to make a map to show the thickness of icesheet.
First try with ListDensityPlot:


You can see it is not working well, we need to limit the plot region inside the boundary of Greenland. Here is the place we can use the point-in-polygon test. In this example, I use inPolyQ2 from the answer by Simon Woods:


Still not right? What's wrong? The trick is to increase MaxPlotPoints  to 100 at least:
ListDensityPlot[data, PlotRange -> All,
 ColorFunction -> (ColorData["Rainbow"][1 - #] &),
 RegionFunction -> (inPolyQ2[greenland_ploygon[[1, 1]], #1, #2] &),
 MaxPlotPoints -> 150, MeshFunctions -> {#3 &}, Mesh -> 10]


It looks like a real map now.

Feb 26, 2014

Fetching data from HTML source


Parsing html to get the data we need can be very frustratingLucky, Mathematica has a powerful hmtl import function, you can import raw html data into several different formats. In my experiences, import html as "XMLObject" is usually the best way to go. 
Here is an example: OSCAR Nominees:
xml = Import["http://oscar.go.com/nominees", "XMLObject"];   
We are interested in the list of nomineed films

body = Cases[xml, XMLElement["div", {"class" -> "nominee-by-film"}, ___], Infinity];
Extract titles:
title = Cases[body, XMLElement["span", {"class" -> "title"}, value_] :> value, Infinity] 
Extract the number of nominees:
nominee =
  Cases[body,
   XMLElement["h1", {"class" -> "numberOfNominations"}, value_] :>
    StringCases[value, x : NumberString :> ToExpression[x]], Infinity] ;
Put these two together:
result = Sort[Transpose[{title, Flatten@nominee}], #1[[2]] > #2[[2]] &]
Let's draw a graph to show the top 10 of the most nomineed films:
oscar = Import["http://www.oscars.org/awards/academyawards/about/awards/images/side_oscar.jpg"];
BarChart[result[[1 ;; 10, 2]],
 ChartLabels -> Placed[Flatten@result[[1 ;; 10, 1]], After],
 BarOrigin -> Left, Background -> LightBlue, ChartElements -> {oscar, {1, 1}},
  Axes -> None, LabelStyle -> {Bold, Darker@Blue, 14}] 

For this particular example, you can also try to get the same information directly from WolframAlpha.

Related post: A discussion on Mathmeatica Stackexchange

Feb 19, 2014

Simple Guide on Geospatial Coordinates Transformation with Mathematica

A few questions on geospatial coordinates transformation have shown up in Mathematica.Stackexchange. Here is a very brief summary.

In US, you probably likely deal with two projection systems: State Plane Coordinates System and UTM.

1. State Plane Coordinates System
In the U.S. State Plane Coordinate System (SPCS), the transverse Mercator projection is used for states that are long in the north-south direction, a Lambert conformal conic projection is used for states that extend in the east-west direction, and the oblique Mercator projection is used for Alaska.

In GeoProjectionData, SPCS83IN01 and SPCS83IN02 represent Indiana Steate Plane east zone and west zone. SPCS83TX01 ~ SPCS83TX05 represent 5 zones from north to south in Texas. Tennessee has only one zone: SPCS83TN00. Here is an online interactive map on SPCS.

There are also SPCS27 series, which are based on NAD27 datum, however, it is quite rare to get the data in the old coordinate system.

One common mistakes is usually caused by the unit: meters vs feet. In Mathematica, the coordinate is calculated in meters, the data you get is probably in feet.

Related posts on stack exchange convert spcs to (lat, lon)convert (lat, lon) to spcs

2. UTM
Universal Transverse Mercator (UTM) coordinate system divides the Earth into sixty zones: UTMZone01 ~ UTMZone60.

In Mathematica 9, there is a problem with UTMZone data:
GeoProjectionData["UTMZone16"]
{"TransverseMercator", {"Centering" -> {0, -87},  "CentralScaleFactor" -> 1, "GridOrigin" -> {0, 0}, "ReferenceModel" -> "WGS84"}}
The scale factor: 0.9996 and the grid origin: {500000,0} shall be specified for coordinate transformation: 
GeoGridPosition[
 GeoPosition[{39.162147, -86.529045}, "WGS84"], {"UTMZone33",
  "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}}]
Related posts on stack exchange convert between (lat, lon) and UTM

This problem is fixed in next version of Mathematica.

Dec 2, 2013

Testing Wolfram Language on a Raspberry Pi emulation

Want to test Wolfram Language without a Raspberry Pi? 

Files you need to download for Windows platform:

  1. QEMU 1.6.0 Binary for Windows: Qemu-1.6.0-windows.zip
  2. Linux Kernel
  3. Latest Raspbian Image 2013-09-25-wheezy-raspbian.zip as this blog is written
  4. The guide on Howto setup Raspberry Pi Emulation with Qemu on Linux or Windows
Main steps:
   1. expand Raspbian image to add more disk space
   2. start the virtual machine to register extra disk space
   3. install Mathematica Language:
       sudo apt-get update && sudo apt-get install wolfram-engine
   
Once you get it running, you can check the performance against a real Raspberry Pi

It takes around 1 hour to get it done. Don't expect much on the performance, have fun!


Updates:


Remote Development Kit doesn't work, the reason is "ssh"


To ssh into the emulated Raspberry Pi, add "-redir tcp:2222::22" to the command options when starting qemu, then "ssh -p 2222 pi@localhost", in Mathematica, it probably connects port 22 by default, it seems no way to specify the port number.

Oct 31, 2011

How to make arrow objects for Google Earth

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

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

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

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

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

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

For the end point:

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

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

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

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

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

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

screenshot