## Weather Maps with Mathematica

22

16

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 :) – Dr. belisarius – 2013-06-03T04:03:44.037

If you want a more detailed map of the USA, this isn't built-in, but there is some data you can use here.

– J. M.'s ennui – 2013-06-03T09:56:26.793

Related, if not duplicate: (528), (2009), (3723)

– Mr.Wizard – 2013-06-03T11:54:02.103

Marginally related, and I may be remembering wrong, but ExampleData[], alone, just by itself, by default, brings up a weather map. – berniethejet – 2013-06-03T14:41:32.447

Is it possible that a Wolfram Blog entry from a few years ago has pointers you can use? The link is http://blog.wolfram.com/2009/01/13/visualizing-weather-patterns-in-mathematica-7/

– dwa – 2013-06-04T07:11:59.147

@dwa, sadly the supporting notebook for that blog entry seems to be gone... (luckily, it can still be obtained from the Wayback Machine.) – J. M.'s ennui – 2013-06-04T09:53:55.123

21

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 – 2013-06-03T08:00:52.180

@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. – J. M.'s ennui – 2013-06-03T08:05:31.417

Dear J. M.,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>>]) – Zviovich – 2013-06-03T14:22:33.610

@Zviovich, hopefully fixed in this version. – J. M.'s ennui – 2013-06-04T01:01:44.040

15

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.