Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
"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. –  J. M. 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

2 Answers 2

up vote 17 down vote accepted

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

share|improve this answer
    
" 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. –  J. M. 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 Graphics`Mesh`InPolygonQ::invpoly: -- Message text not found -- (Graphics`Mesh`InPolygonQ[{{{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. –  J. M. Jun 4 '13 at 1:01
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.