23
$\begingroup$

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

enter image description here

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

enter image description here

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

Improve this question
$\endgroup$
6
  • $\begingroup$ "Pato Criollo" ... ha! Nice username :) $\endgroup$ Commented Jun 3, 2013 at 4:03
  • $\begingroup$ If you want a more detailed map of the USA, this isn't built-in, but there is some data you can use here. $\endgroup$ Commented Jun 3, 2013 at 9:56
  • $\begingroup$ Related, if not duplicate: (528), (2009), (3723) $\endgroup$ Commented Jun 3, 2013 at 11:54
  • $\begingroup$ Marginally related, and I may be remembering wrong, but ExampleData[], alone, just by itself, by default, brings up a weather map. $\endgroup$ Commented Jun 3, 2013 at 14:41
  • $\begingroup$ 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/… $\endgroup$ Commented Jun 4, 2013 at 7:11

2 Answers 2

Reset to default
21
$\begingroup$

One key function you might need is the (undocumented) function Graphics`Mesh`InPolygonQ[], 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, Graphics`Mesh`InPolygonQ[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"]

heat map for USA

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 @@ (Graphics`Mesh`InPolygonQ[#, {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"]

USA temperature map with finer sampling

The color function can be easily changed, of course:

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

USA temperature map with finer sampling, different color scheme

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

weatherMap["WesternEurope", "Temperature"]

temperature map for Western Europe

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]

Malaysian temperature map, jet color scheme

(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.)

Improve this answer
$\endgroup$
4
  • $\begingroup$ " 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? $\endgroup$ Commented Jun 3, 2013 at 8:00
  • $\begingroup$ @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. $\endgroup$ Commented Jun 3, 2013 at 8:05
  • $\begingroup$ 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>>]) $\endgroup$ Commented Jun 3, 2013 at 14:22
  • $\begingroup$ @Zviovich, hopefully fixed in this version. $\endgroup$ Commented Jun 4, 2013 at 1:01
15
$\begingroup$

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

result

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.

Improve this answer
$\endgroup$
0

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.