# Weather Maps with Mathematica

I've been trying to set up a weather map with Mathematica similar to what NOAA publishes in their website.

So far I've been successful to bring the data over and extend it across the whole picture frame.

weatherMap[region_, property_] :=
Module[{polygon, coord, minX, maxX, minY, maxY, data, deltaX, deltaY},

polygon = CountryData[region, "Polygon"];
coord = Cases[CountryData[region, "Coordinates"], {_, _}, Infinity];
minX = Floor@Min[coord[[All, 2]]];
maxX = Ceiling@Max[coord[[All, 2]]];
minY = Floor@Min[coord[[All, 1]]];
maxY = Ceiling@Max[coord[[All, 1]]];
deltaX = (maxX - minX)/25;
deltaY = (maxY - minY)/25;
data = DeleteCases[
Flatten[Table[{x, y, WeatherData[{y, x}, property]}, {x, minX,
maxX, deltaX}, {y, minY, maxY, deltaY}], 1], {x_,
y_, _Missing}];
Show[Graphics[{Opacity[0], polygon}],
ListContourPlot[data, Contours -> 10,
ColorFunction -> "Temperature", InterpolationOrder -> 3,
PlotLegends -> Placed[Automatic, Above],
PerformanceGoal -> "Quality"],
Graphics[{EdgeForm[Black], FaceForm[Opacity[0.01], GrayLevel[.1]],
polygon}]]]
weatherMap["USA", "Temperature"]

What would be the best way to limit the contour temperature graphic to within the frame of the given polygon?

-
"Pato Criollo" ... ha! Nice username :) –  belisarius Jun 3 '13 at 4:03
If you want a more detailed map of the USA, this isn't built-in, but there is some data you can use here. –  Ｊ. Ｍ. Jun 3 '13 at 9:56
Related, if not duplicate: (528), (2009), (3723) –  Mr.Wizard Jun 3 '13 at 11:54
Marginally related, and I may be remembering wrong, but ExampleData[], alone, just by itself, by default, brings up a weather map. –  berniethejet Jun 3 '13 at 14:41
Is it possible that a Wolfram Blog entry from a few years ago has pointers you can use? The link is blog.wolfram.com/2009/01/13/… –  dwa Jun 4 '13 at 7:11
show 1 more comment

One key function you might need is the (undocumented) function GraphicsMeshInPolygonQ[], which tests if a point is inside a given polygon. With it, and a few other tweaks, here's my version of weatherMap[]:

weatherMap[region_String, property_String, res_Integer: 25, opts___] :=
Module[{fmin, cmax, coords, pts, minLong, maxLong, minLat, maxLat,
deltaLong, deltaLat, data},
fmin = Composition[Floor, Min]; cmax = Composition[Ceiling, Max];
coords = CountryData[region, "Coordinates"] /. v_?VectorQ :> Reverse[v];
pts = Flatten[coords, 1];
{minLong, maxLong} = Through[{fmin, cmax}[pts[[All, 1]]]];
{minLat, maxLat} = Through[{fmin, cmax}[pts[[All, 2]]]];
{deltaLong, deltaLat} = {maxLong - minLong, maxLat - minLat}/res;
data = Table[WeatherData[{lat, long}, property],
{long, minLong, maxLong, deltaLong}, {lat, minLat, maxLat, deltaLat}];
ListContourPlot[data, AspectRatio -> Automatic,
DataRange -> {{minLong, maxLong}, {minLat, maxLat}},
Epilog -> {Directive[EdgeForm[Black], FaceForm[None]], Polygon[coords]},
opts, ColorFunction -> "ThermometerColors", Contours -> 10,
InterpolationOrder -> 3, PerformanceGoal -> "Quality",
RegionFunction -> (Or @@ Map[Function[polys,
GraphicsMeshInPolygonQ[polys, {#1, #2}]], coords] &)]]

(Not having Mathematica 9, I had to omit the option PlotLegends; you can add it back yourself.)

Test:

weatherMap["USA", "Temperature"]

I should note that I thought it best not to delete Missing[] entries, since ListContourPlot[] is certainly able to cope with them (note the holes in the picture), and I feel it is dishonest to pretend that data exists in regions where there is actually none.

If you want your weather maps to be a bit more accurate, and you don't mind a bit of a wait (WeatherData[] is, after all, rather slow), you can use a DensityPlot[] (similar to what Mark suggested), and let the internal adaptive sampling algorithms do the sampling for you, instead of having to sample at equispaced longitudes/latitudes, and potentially missing data. With this in mind, here is a shorter (but slower) implementation of weatherMap[]:

weatherMap[region_String, property_String, opts___] :=
Module[{coords, minLong, maxLong, minLat, maxLat},

coords = CountryData[region, "Coordinates"] /. v_?VectorQ :> Reverse[v];
If[VectorQ[CountryData[region, "Countries"], TrueQ], coords = Flatten[coords, 1]];
{{minLong, maxLong}, {minLat, maxLat}} = PlotRange[Graphics[Polygon /@ coords]];

DensityPlot[WeatherData[{lat, long}, property],
{long, minLong, maxLong}, {lat, minLat, maxLat},
AspectRatio -> Automatic, MeshFunctions -> {#3 &}, opts,
BoundaryStyle -> Directive[Thickness[Medium], Black],
ColorFunction -> "ThermometerColors", MaxRecursion -> 1,
Mesh -> 10, PerformanceGoal -> "Quality",
RegionFunction -> Function[{x, y, z},
Or @@ (GraphicsMeshInPolygonQ[#, {x, y}] & /@ coords)]]]

Some new features you might notice include the (undocumented) use of PlotRange[] to reckon out appropriate bounds for the data to be plotted, and the use of MeshFunctions -> {#3 &} so that one could still see the contours where property is constant.

Here again is the temperature map for the United States:

weatherMap["USA", "Temperature"]

The color function can be easily changed, of course:

weatherMap["USA", "Temperature", ColorFunction -> "LightTemperatureMap"]

Here's a demonstration of weatherMap[] on a collection of countries, showing its handling of noncontiguous land masses:

weatherMap["WesternEurope", "Temperature"]

Finally, to demonstrate the use of a "rainbow" map, here's a temperature map of Malaysia:

(* "jet" colormap *)
jet[u_?NumericQ] :=
Blend[{{0, RGBColor[0, 0, 9/16]}, {1/9, Blue}, {23/63, Cyan}, {13/21, Yellow},
{47/63, Orange}, {55/63, Red}, {1, RGBColor[1/2, 0, 0]}}, u] /; 0 <= u <= 1

weatherMap["Malaysia", "Temperature", ColorFunction -> jet]

(On the other hand, there are a number of reasons to avoid rainbow-like color maps for applications like these; see e.g. this and this.)

-
" it is dishonest to pretend that data exists in regions where there is actually none." Isn't that just what we'd usually call interpolation? –  Sjoerd C. de Vries Jun 3 '13 at 8:00
@Sjoerd, sure, but if there are no measurements within a large swath of the region being considered, interpolation isn't exactly honest; interpolation might say a patch the size of Detroit is currently sunny if the surrounding counties are sunny, even though there is in fact a nice spell of sleet currently happening. –  Ｊ. Ｍ. Jun 3 '13 at 8:05
Dear 0x4A4D,There are some issues with the code when trying to run it against a group of countries, such as weatherMap["WesternEurope","Temperature"]. Any ideas on how to solve it? Error message GraphicsMeshInPolygonQ::invpoly: -- Message text not found -- (GraphicsMeshInPolygonQ[{{{1.72194,42.5012},{1.70363,42.4898},{1.67035,42.5041‌​},{1.65746,42.4687},{1.59763,42.4676},{1.58084,42.45},<<16>>,{1.5784,42.6461},{1.‌​60504,42.6231},{1.74013,42.6129},{1.73606,42.5875},{1.79514,42.5821},{1.73871,42.‌​5655}}},<<1>>]) –  PatoCriollo Jun 3 '13 at 14:22
@Pato, hopefully fixed in this version. –  Ｊ. Ｍ. Jun 4 '13 at 1:01

If you want it to look like your NOAA example, that's not a ContourPlot[]. It's a DensityPlot[] using a modified Hue[] function for the ColorFunction. E.g. something like:

Show[ListDensityPlot[{#4, #3, #5} & @@@ Select[maxt, #[[5]] < 9999 &],
ColorFunctionScaling -> False,
ColorFunction -> Function[{x}, Hue[Clip[(100 - x)/60, {0, 1}]]],
AspectRatio -> Automatic],
Graphics[{EdgeForm[Black], Opacity[0],
CountryData["UnitedStates", "Polygon"]}]]

See the other answer for a polygon solution.

Update:

@0x4A4D asked how to get maxt. I downloaded the data from NOAA, ds.maxt.bin, which is in GRIB2 format. Unfortunately Mathematica 9's Import[] knows about GRIB, but not GRIB2. So I also downloaded, compiled, and ran code to translate GRIB2 into CSV. I picked one of the three data sets generated, imported that, and assigned it to maxt.

-