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.

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

Control Google Earth from Mathematica

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

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

SC[{lat_,lon_}]:=r {Cos[lon \[Degree]] Cos[lat \[Degree]],Sin[lon \[Degree]] Cos[lat \[Degree]],Sin[lat \[Degree]]};
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.

Create the map with a image background

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]

Tip: create 3D Map

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.

 

Nationality Name

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.

Fun with Flags: who have star in their flags?

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.

Leader of the world? again

Use the flags at the vertices of the graph.
VertexRenderingFunction -> (Inset[Show[CountryData[#2, "Flag"], ImageSize -> 30], #] &)

See the original.

Leader of the World?

The graphplot of all the countries with their major export partners.

GraphPlot[DeleteCases[Flatten[Thread[# -> CountryData[#, "ExportPartners"] /. _Missing :> {}] & /@ CountryData[]], _ -> {}], VertexLabeling -> True, MultiedgeStyle -> True, DirectedEdges -> True, Method -> HighDimensionalEmbedding]

Playing with numbers: Oil Consumption

Ok, everyone knows the fact:

oil consumption in barrels per day

(It's probably the data of year 2005)

G8 vs China