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 am trying to use ListVectorPlot in combination with WeatherData and CountryData to show wind velocity vectors over a given country. The problem I am having stems from "NotAvailable" and "NotApplicable" data from a function in which the wind's velocity vector is constructed. The function is

WindFoo[x_, y_] := -WeatherData[{y, x}, "WindSpeed"] {1/Cos[y °], 
   1} Through[{Sin, Cos}[WeatherData[{y, x}, "WindDirection"] °]]

Note: I am not able to display the degree symbol, so it is written out in the above code.

How can the code above be modified to exclude the cases of "NotAvailable" and "NotApplicable" data to work in manner similar to the wind velocity vectors shown over Australia in the following tutorial: http://blog.wolfram.com/2009/01/13/visualizing-weather-patterns-in-mathematica-7/

share|improve this question
    
what version of Mathematica are you using? –  ubpdqn Jul 27 at 3:02
1  
@ubpdqn version 10 –  sunspots Jul 27 at 3:11
    
GeoGraphics in version 10 can be used for map but requires some work to get the vector plot as you cannot combine with Show. However, you can use CountryData["Australia","Schematic Polygon"]. With respect to the removal of missing data, there are two issues. In v 10 the above produces things with units so need to strip QuantityMagnitude and probably easier to do Cases where data all numeric. –  ubpdqn Jul 27 at 3:16

3 Answers 3

up vote 6 down vote accepted

This uses some version 10 functions:

Getting the data (takes time as I could have set this up better):

dat = Table[{i, j, 
    WeatherData[{j, i}, #] & /@ {"WindSpeed", "WindDirection"}}, {i, 
    113, 153, 4}, {j, -43, -11, 4}];

Processing:

data = Cases[{{#1, #2}, -QuantityMagnitude[#3[[1]]] {1/Cos[#2 Degree], 
        1} Through[{Sin, Cos}[#3[[2]]]]} & @@@ 
    Flatten[dat, 1], {_, {_?NumericQ, _?NumericQ}}];

Visualizing:

ls = ListStreamPlot[data];
arrows = Cases[ls, Arrow[_], Infinity];
GeoGraphics[{Purple, arrows}]

enter image description here

Please look at documentation as must reverse coordinates if you choose to use GeoPosition. Obviously can modify to suit stylistic preferences and other needs.

UPDATE

Just to 'check' this with another source (noting I am not sure what WeatherData uses as default, I just assumed it used the most recent data): From http://www.windfinder.com/weather-maps/forecast/australia#4/-24.69/\ 135.79 enter image description here

or a little before:

enter image description here

These seem somewhat similar to the Mathematica generated map.

Exploiting the hard work of Gleno and appropriately adapting code from the documentation of WindVectorData I believe yields the "correct" wind stream plot consistent with the check above (the issue is transforming the wind velocity components):

windVelocity[
  coords_] := {Reverse@
   coords, -Reverse@QuantityMagnitude[WindVectorData[coords]] {1/
     Cos[coords[[1]] \[Degree]], 1}}

wvd[country_String, date_, n_] := 
  Module[{box, latmin, latmax, longmin, longmax, pts, sample, stream, 
    latstep, longstep, ls, arrows}, 
   box = {Min@#, Max@#} & /@ 
     Transpose@
      Cases[CountryData[country, "Polygon"], {_Real, _}, Infinity];
   {{latmin, latmax}, {longmin, longmax}} = box + {{-2, 2}, {-2, 2}};
   latstep = (latmax - latmin)/(n - 1);
   longstep = (longmax - longmin)/(n - 1);
   pts = Flatten[
     Table[{lat, long}, {lat, latmin, latmax, latstep}, {long, 
       longmin, longmax, longstep}], 1];
   sample = windVelocity /@ pts;
   stream = Cases[sample, u_ /; FreeQ[u, "NotAvailable"]];
   ls = ListStreamPlot[stream];
   arrows = Cases[ls, Arrow[_], Infinity];
   GeoGraphics[arrows, GeoBackground -> "ReliefMap"]
   ];

This leads to:

gr = wvd["Australia", Today, 10]

enter image description here

I guess the answer my friend is blowing in the wind...

share|improve this answer
    
The new maps in Mathematica 10 look really cool. I think at least one of us botched something, because we get quite stark disagreement on wind direction. I'll check tomorrow where the wind actually blows. ;) –  Gleno Jul 27 at 4:26
    
@Gleno its probably me. I am just checking. Very much like yuor answer. Gets my vote. –  ubpdqn Jul 27 at 4:30
    
@Gleno see my update...I think I have resolve this issue –  ubpdqn Jul 27 at 6:22
    
Ah, well done sir! –  Gleno Jul 27 at 16:06
1  
I don't think I ever congratulated you on reaching 10K "reputation." Nice job! It's been fun watching the quality of your answers improve over time. I like to think I played some role in that. :-) I know how long it takes to write 425 answers; thanks for giving back to the community! –  Mr.Wizard Jul 27 at 23:29

NB I think there's a problem with this code, please see this answer by ubpdqn for more information. I'll update these snippets when I get the chance.

Yours is a question with many possible interpretations. I've chosen the interpretation that was most fun for me to play with, so...

ang = 20; (* divide the world into chunks of this size *)
pts = Flatten[Table[{lat, long}, {lat, -80, 80, ang}, {long, -180, 180, ang}], 1];
sample = WindVectorData[Replace[pts, c_ :> {c, Today}, 1]]; (* this takes forever *)

stream = Cases[Transpose@{Reverse/@pts, sample}, {_, Except[_Missing]}] /. 
    Quantity[v_, _] :> v /. Interval[lst_] :> Mean /@ lst;

(* This function crashes kernel very often, but in interpolation, I assume *)
streamplot = ListStreamPlot[stream, StreamColorFunction -> "SolarColors", 
   InterpolationOrder -> None, StreamPoints -> Fine];

mpoly[poly_] := poly /. GeoPosition[pts_] :> pts 
                     /. {lat_Real, long_} :> {long, lat} 
                     /. Polygon[data_] :> Style[Polygon[data], Antialiasing -> True];

world = Graphics[{RGBColor[0.38, 0.77, 1.]}~Join~
    Table[mpoly@CountryData[c, "SchematicPolygon"], {c, CountryData[]}]];
g = Show[world, streamplot, ImageSize -> 1400]

If you have Mathematica 10, the graph g should be:

Worlwinds

Puting it together, you can have a function like

PlotWindBox[country_String, date_, n_] := 
  Module[{box, latmin, latmax, longmin, longmax, pts, sample, stream, 
    streamplot, mpoly, target, world, latstep, longstep},

   box = {Min@#, Max@#} & /@ 
     Transpose@
      Cases[CountryData[country, "Polygon"], {_Real, _}, Infinity];
   {{latmin, latmax}, {longmin, longmax}} = box + {{-2, 2}, {-2, 2}};


   latstep = (latmax - latmin)/(n - 1);
   longstep = (longmax - longmin)/(n - 1);
   pts = Flatten[
     Table[{lat, long}, {lat, latmin, latmax, latstep}, {long, 
       longmin, longmax, longstep}], 1];

   sample = WindVectorData[Replace[pts, c_ :> {c, date}, 1]];
   stream = 
    Cases[Transpose@{Reverse /@ pts, sample}, {_, 
        Except[_Missing]}] /. Quantity[v_, _] :> v /. 
     Interval[lst_] :> Mean /@ lst;
   streamplot = 
    ListStreamPlot[stream, StreamColorFunction -> "SolarColors", 
     StreamPoints -> Fine, InterpolationOrder -> None];
   mpoly[poly_] := 
    poly /. 
       GeoPosition[pts_] :> pts /. {lat_Real, long_} :> {long, 
        lat} /. Polygon[data_] :> 
      Style[Polygon[data], Antialiasing -> True];
   target = 
    Graphics[{RGBColor[0.38, 0.77, 1.], 
      mpoly@CountryData[country, "Polygon"]}];
   world = 
    Graphics[{GrayLevel[0.84]}~Join~
      Table[mpoly@CountryData[c, "Polygon"], {c, 
        Cases[CountryData[], Except[Entity["Country", country]]]}]];
   Show[world, target, streamplot, ImageSize -> 1400, 
    PlotRange -> {{longmin, longmax}, {latmin, latmax}}]
   ];

Call it with PlotWindBox["Australia", Today, 10] (10 means that there are 10^2 sample points, try with fewer sample points if you are impatient like me) and obtain:

Australia Weather

share|improve this answer

I feel a bit foolish posting such a simple answer after the two complete examples above, but I think it bears pointing out that version 10 includes a function specifically for removing Missing values:

DeleteMissing:

enter image description here

share|improve this answer
    
version 10 keeps coming up with surprises some nice, like this one...very handy so thanks...many I am still adjusting to –  ubpdqn Jul 27 at 11:35

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.