2
$\begingroup$

I am currently working on a project where I am trying to plot latitude, longitude, and maximum rainfall of hurricane Harvey on a 3D plot. I have gathered all the data inside Mathematica using station names and the WeatherData function. I am trying to use ListInterpolation to smooth the graph I am obtaining from data, but the problem is that my data contains latitude and longitude sets as angles in degrees. This is the code I use to obtain my data:

coordinates = CityData[{"Dallas", "Texas", "UnitedStates"}, "Coordinates"];
weatherstations = WeatherData[{coordinates, 400}];
Lat = Table[WeatherData[weatherstations[[i]], "Latitude"], {i, 400}];
Lon = Table[WeatherData[weatherstations[[j]], "Longitude"], {j, 400}];
dRain = 
  Table[
    WeatherData[
      weatherstations[[k]], "TotalPrecipitation", 
      {{2017, 8, 17}, {2017, 9, 3}, "Day"}], 
    {k, 400}];
MaxRain = Table[DeleteCases[Max[dRain[[k]]], _Missing], {k, 400}];
Rdata = DeleteCases[Transpose[{Lat, Lon, MaxRain}], {_, _, _Missing}];

enter image description here

The data can now be plotted using 3DListPlot; however, I cant use ListInterpolation[data] because the data contains degree symbols. Does anyone know how do I get around this problem? Do you know what I am trying to say?

$\endgroup$
  • 1
    $\begingroup$ ListInterpolation is almost certainly the wrong function to use; you'll probably want to try Interpolation instead. Could you please update your answer with actual code rather than an image to make it easier to copy? That makes it a lot easier to help out. $\endgroup$ – Sjoerd Smit Nov 17 '17 at 14:31
  • $\begingroup$ I tried to use Interpolation and I came up with an error as well. I can edit and post the code right now for you. Thank you. $\endgroup$ – hwhorf Nov 17 '17 at 14:35
  • $\begingroup$ Interpolation asks for data in a different format than the one you seem to be using. I'm also not sure if it will accept Quantitys, so you may have to convert the coordinates to regular numbers. $\endgroup$ – Sjoerd Smit Nov 17 '17 at 14:36
  • $\begingroup$ Converting the coordinates to regular numbers? Like degrees to radians? Or is there a way to drop the degree symbol. Btw, to run the code I don't recommend using 400 data points like what I did. $\endgroup$ – hwhorf Nov 17 '17 at 14:45
3
$\begingroup$

My answer also relies on retrieving the property "Coordinates" to get the coordinates as real numbers, but it also shows how the code can be written without using "Table" over and over again.

dallas = CityData[{"Dallas", "Texas", "UnitedStates"}, "Coordinates"];

With[{n = 400}, weatherStations = WeatherData[{dallas, n}]];
coordinates = Through[weatherStations["Coordinates"]];
rain =
  WeatherData[#, "TotalPrecipitation", {{2017, 8, 17}, {2017, 9, 3}, "Day"}] & /@ 
    weatherStations;
maxRain =
  (Max /@ DeleteCases[Through[rain["Values"]], Missing[_], {-2}, Heads -> True]) /. 
    {"Values" -> 0., -∞ -> 0.};
data = MapThread[{Sequence @@ #1, #2} &, {coordinates, maxRain}];

Now we have data that is all cleaned up and plots nicely. The red sphere floats above the location of Dallas as retrieved in the code.

Show[
  ListPlot3D[data, PlotRange -> All, BoxRatios -> {1, 1, 1/2}],
  Graphics3D[{Red, Ellipsoid[{32.7942, -96.7655, 2.}, .015 {16, 18, 80}]}],
  ImageSize -> 600]

plot

The plot appears to show that the heaviest rainfall was concentrated along a relative short line segment located south-by-southeast of the city center.

Update

This update addresses an issue raised by the OP in a comment to this answer.

The maximum wind speed at each weather station can be retrieved in pretty much the same way as the rain data, and it can be used to colorize the rain data. Doing so does add a fair amount more work to making the plot.

Here is how I did it. Note that I have changed the rain data dependent variable from amount of rain on the day of it rained hardest to the total accumulation of rain over the interval under consideration. It doesn't have much effect on the shape of the plot, but it does increased the magnitude of the quantities plotted. You can ignore this change if you wish—just substituted data from the 1st part of this answer for rainData.

dallas = CityData[{"Dallas", "Texas", "UnitedStates"}, "Coordinates"];
With[{n = 400}, weatherStations = WeatherData[{dallas, n}]];
coordinates = Through[weatherStations["Coordinates"]];
dailyRain =
  WeatherData[#, "TotalPrecipitation", {{2017, 8, 17}, {2017, 9, 3}, "Day"}] & /@ 
    weatherStations;
totalRain =
  Total[
    (DeleteCases[Through[dailyRain["Values"]], Missing[_], {-2}, Heads -> True]) /. 
      {"Values" -> {0}, {} -> {0}}, 
    {2}];
rainData = MapThread[{Sequence @@ #1, #2} &, {coordinates, totalRain}];
dailyMaxWinds =
  WeatherData[#, "MaxWindSpeed", {{2017, 8, 17}, {2017, 9, 3}, "Day"}] & /@ 
    weatherStations;
maxWinds = 
  (Max /@ 
     DeleteCases[Through[dailyMaxWinds["Values"]], Missing[_], {-2}, Heads -> True]) 
   /. {"Values" -> 0., -∞ -> 0.};
(* 
   the following associates each coordinate with a maxWinds value scaled to lie 
   in the unit interval 
*)
windData = MapThread[{#1, #2} &, {coordinates, maxWind/Max[maxWinds]}];
windF = Quiet @ Interpolation[windData];

Legended[
  Show[
    ListPlot3D[rainData,
      (* this is where windF is used to colorize the plot *)
      ColorFunction -> (ColorData["ThermometerColors"][windF[#1, #2]] &),
      ColorFunctionScaling -> False,
      PlotRange -> All,
      BoxRatios -> {1, 1, 1/2}],
    Graphics3D[{Red, Ellipsoid[{32.7942, -96.7655, 3.}, .015 {16, 18, 2 74}]}],
    ImageSize -> 600],
  BarLegend[{"ThermometerColors", {0., Max[maxWinds]}}, 5]]

plot

$\endgroup$
  • $\begingroup$ A, nice application of Through. I always kind of forget when to use that one. $\endgroup$ – Sjoerd Smit Nov 17 '17 at 16:57
  • $\begingroup$ This is awesome! However, the graph looks a bit too...orange. Lets say I want to color code this graph based on the values of "WindSpeed". This means I want higher values of wind speed red and lower values orange. You could gather your data points for wind speed the same as how you gathered them for "TotalPercipitation"( just by using WeatherData) but how do you add this extra dimension to the graph? $\endgroup$ – hwhorf Nov 17 '17 at 18:31
  • $\begingroup$ @HarrisonWhorf. I have updated my answer to address the matter of coloring the surface to show maximum wind velocity. $\endgroup$ – m_goldberg Nov 18 '17 at 8:22
3
$\begingroup$

You're pretty close, but the data format you're constructing will not work for Interpolate. This should get you where you need to be:

n = 20;
coordinates = 
  CityData[{"Dallas", "Texas", "UnitedStates"}, "Coordinates"];
weatherstations = WeatherData[{coordinates, n}];
coordinates = 
  Table[WeatherData[weatherstations[[i]], "Coordinates"], {i, n}];
dRain = Table[
   WeatherData[weatherstations[[k]], 
    "TotalPrecipitation", {{2017, 8, 17}, {2017, 9, 3}, "Day"}], {k, 
    n}];
MaxRain = Table[DeleteCases[_Missing]@Max[dRain[[k]]], {k, n}];
Rdata = DeleteCases[Transpose[{coordinates, MaxRain}], {_, _Missing}];

Note that the coordinates are now regular numbers instead of Quantity objects (you can convert quantities to numbers with QuantityMagnitude). You can now generate an interpolation function:

int = Interpolation[Rdata]

Note the message that warns you that the interpolation will be 1st order because of the irregular spacing of the grid points. If you want to have a smoother interpolation, you can try to see if Predict works on your data:

pred = Predict[Rule @@@ Rdata]

Note that these 2 work differently if you want to make a prediction. E.g., if you want to predict for the coordinates {32, -97}, you do:

int[32, -97]

or

pred[{32, -97}]
$\endgroup$
  • $\begingroup$ Thank you. Using "Coordinates" will give you lat, lon as actual numbers. Predict still doesn't work on my data however. $\endgroup$ – hwhorf Nov 17 '17 at 15:28
  • $\begingroup$ Because my data is irregularly spaced, can these functions still be used?Perhaps a better way to visualize this data is to plot these point on a map of Texas and color code them based on the value of height. Do you know how do that? The function GeoGraphics[{GeoBoundsRegion[{{Min[Lat], Max[Lat]}, {Min[Lon], Max[Lon]}}]}] will display the bounds of my data. $\endgroup$ – hwhorf Nov 17 '17 at 15:34
  • $\begingroup$ Well, you don't need to worry too much about the irregular spacing. It's just that the interpolation will be 1st order, so it'll be a bit blocky. And Predict doesn't work for you? What version of Mathematica are you on? I'm on 11.2 and if I copy the code above verbatim into a notebook, everything works. If you want to use GeoGraphics, I recommend you look at GeoBubbleChart (if your version has it). Finally, you can also visialise the data directly with something like ListPlot3D[Flatten /@ Rdata] $\endgroup$ – Sjoerd Smit Nov 17 '17 at 15:42
  • $\begingroup$ I am on version 11.0. I am going to update to version 11.2 right now. Just thinking of it, could I use linspace (or another function for that matter) to generate an evenly spaced grid that is the bounds of my longitude and latitude. MATLAB's interpolation function allows you to input your irregularly spaced data and the evenly spaced meshgrid to "smooth" out the graph. Is there a way to do this in Mathematica? $\endgroup$ – hwhorf Nov 17 '17 at 16:00
  • $\begingroup$ Well, you could do something like sampleGrid = Flatten[CoordinateBoundsArray[{{minX, maxX}, {minY, maxY}}, Into[20]], 1] for some suitable values for minX etc. Then you can resample the function values on that grid with int @@ Transpose[sampleGrid] or pred[sampleGrid]. $\endgroup$ – Sjoerd Smit Nov 17 '17 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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