Message Boards Message Boards

Climate Change in the United States

The U.S. National Oceanic and Atmospheric Administration produces large datasets containing long-term time series of weather data for various geographical regions of the United States. I decided to use the most recent file of monthly average temperatures ("climdiv-tmpcst-v1 .0.0-20181204") to produce a map of climate change by State. The code works, but as there are often many ways to do the same thing in the Wolfram Language, I would be very interested to know if anyone can suggest faster, more compact or more readable methods of doing what I am doing here.

I begin by reading in the file,which consists of 13 space-delimited fields. The first encodes the region and year, and the next 12 are monthly average temperatures in degrees Fahrenheit. Additional information about the data is provided in ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv/state-readme.txt.

tdata = SemanticImport[
   "ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv/climdiv-tmpcst-v1.0.\
0-20181204", {"String", "Number", "Number", "Number", "Number", 
    "Number", "Number", "Number", "Number", "Number", "Number", 
    "Number", "Number"}];

The file has no headers, so I add some:

tdata = tdata[
   All, <|"Code" -> 1, "Jan" -> 2, "Feb" -> 3, "Mar" -> 4, 
    "Apr" -> 5, "May" -> 6, "Jun" -> 7, "Jul" -> 8, "Aug" -> 9, 
    "Sep" -> 10, "Oct" -> 11, "Nov" -> 12, "Dec" -> 13|>];

Next I add two columns to the dataset: "Region", which corresponds to the first three characters of the Code field, and "Year", which corresponds the last four digits of the Code:

tdata = tdata[
   All, #~Join~<|"Region" -> ToExpression[StringTake[#Code, 3]], 
      "Year" -> ToExpression[StringTake[#Code, -4]]|> &];

I then add "Tavg", the average annual temperature, which is a weighted average of the monthly temperatures. For the sake of accuracy I also account for leap years:

tdata = tdata[
   All, #~Join~<|
      "Tavg" -> (31*#Jan + 31*#Mar + 30*#Apr + 31*#May + 30*#Jun + 
          31*#Jul + 31*#Aug + 30*#Sep + 31*#Oct + 30*#Nov + 
          31 + #Dec + If[LeapYearQ[#Year], 29, 28]*#Feb)/
        If[LeapYearQ[#Year], 366, 365]|> &];

Region codes 1 through 48 correspond to the contiguous States. I am only interested in these, and in the years prior to 2018 (since the data is not yet complete for last year):

tdata = tdata[Select[#Region <= 48 && #Year < 2018 &]];

Next, I extract just the region code, the year and the average temperature:

tdata = tdata[All, {"Region", "Year", "Tavg"}];

This is a list of the 48 contiguous states:

contiguous = {"Alabama", "Arizona", "Arkansas", "California", 
   "Colorado", "Connecticut", "Delaware", "Florida", "Georgia", 
   "Idaho", "Illinois", "Indiana", "Iowa", "Kansas", "Kentucky", 
   "Louisiana", "Maine", "Maryland", "Massachusetts", "Michigan", 
   "Minnesota", "Mississippi", "Missouri", "Montana", "Nebraska", 
   "Nevada", "NewHampshire", "NewJersey", "NewMexico", "NewYork", 
   "NorthCarolina", "NorthDakota", "Ohio", "Oklahoma", "Oregon", 
   "Pennsylvania", "RhodeIsland", "SouthCarolina", "SouthDakota", 
   "Tennessee", "Texas", "Utah", "Vermont", "Virginia", "Washington", 
   "WestVirginia", "Wisconsin", "Wyoming"};

Given this list, we can add a new column containing Entity values for each of the States:

tdata = tdata[
   All, #~Join~<|
      "StateName" -> 
       Entity["AdministrativeDivision", {contiguous[[#Region]], 
         "UnitedStates"}]|> &];

There are of course any number of ways to characterize climate change. First, since annual average temperatures are somewhat noisy we need some type of smoothing. Here I chose a moving average with a window of ten years and plotted the change in the smoothed temperature over the past 100 years:

f[s_List] := 
 Module[{ma}, ma = MovingAverage[s, 10]; (Last[ma] - ma[[-100]])]

s = tdata[GroupBy["StateName"], f, "Tavg"];

GeoRegionValuePlot[s]

Change in average annual temperature over the past century in degrees Fahrenheit

The figure shows that over the past century, the southern states experienced the least warming, with more warming as one proceeds outward to the north, west and east. This is generally in line with other maps I have seen.

Questions: How do I put a title on this map? How do I label the legend?

Attachments:
climatechange.nb
POSTED BY: John Shonder
12 Replies
Sort By:

Does this help point you in the right direction?

legend = SwatchLegend[
   Automatic, {"0\[Degree]C", "1\[Degree]C", "2\[Degree]C", 
    "3\[Degree]C"}, 
   LegendMarkers -> 
    Graphics[{EdgeForm[Black], Opacity[1], Rectangle[]}], 
   LegendLabel -> "Temp \[Degree]C", 
   LegendFunction -> (Framed[#, RoundingRadius -> 5] &), 
   LegendMargins -> 5];
GeoRegionValuePlot[s, Frame -> True, FrameTicks -> None, 
 FrameLabel -> {"Bottom Title", "Left Title", "Top Title", 
   "Right Title"}, PlotLegends -> Placed[legend, Right]]

Title Climate Map
POSTED BY: Tim Laska

This representation implies that the climate obeys state boundaries, which is nonsense.

In my experience maps like this are quite common, and most people know how to interpret them. The one below, for example. Despite the pronounced differences between states, I don't think anyone would conclude that tornadoes obey state boundaries.

Tornadoes per Square Mile by State

As for the rise in annual temperature, given data at a smaller scale one would not see such sharp boundaries between the states. For example, here is a map of temperature rise by Climate Division over a similar period. It looks to me like it is telling the same story as the map I made: marginal warming in the Deep South, and more warming in the East and the West. At this scale there are fewer sharp differences at state boundaries, but perhaps you will now conclude that climate obeys the boundaries of NOAA's climate divisions.

Temperature Rise by US Climate Division
POSTED BY: John Shonder
Rohit Namjoshi
Rohit Namjoshi
Posted 6 years ago

Hi John,

Nice analysis. I had to make a couple of changes to get it to work. Looks like the NOAA just has the most recent data on their FTP server. The latest data as of today is at ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv/climdiv-tmpcst-v1.0.0-20190408. I also had to change all of the Number in the SemanticImport to Real to avoid an odd error message:

MapThread::mptd: Object
Interpreter[Number][{0.00,-0.10,-0.20,0.20,-0.30,0.30,-0.40,0.40,-0.50,-0.60,0.60,-0.70,0.70,-0.80,0.80,-0.90,0.90,-1.00,1.00,10.00,-10.10,10.10,10.20,10.30,10.40,10.50,10.60,10.70,10.80,-10.90,10.90,1.10,11.00,11.10,11.20,11.30,11.40,11.50,11.60,11.70,-11.80,11.80,-11.90,11.90,-1.20,1.20,12.00,-12.10,12.10,12.20,<<731>>}] at position {2, 2} in
MapThread[Rule,{{0.00,-0.10,-0.20,0.20,-0.30,0.30,-0.40,0.40,-0.50,-0.60,0.60,-0.70,0.70,-0.80,0.80,-0.90,0.90,-1.00,1.00,10.00,-10.10,10.10,10.20,10.30,10.40,10.50,10.60,10.70,10.80,-10.90,10.90,1.10,11.00,11.10,11.20,11.30,11.40,11.50,11.60,11.70,-11.80,11.80,-11.90,11.90,-1.20,1.20,12.00,-12.10,12.10,12.20,<<731>>},Interpreter[Number][{<<1>>}]}]

has only 0 of required 1 dimensions.

Some observations

There is no need for ~Join~, associations implicitly join so this

tdata = tdata[
All, #~Join~<|"Region" -> ToExpression[StringTake[#Code, 3]], 
"Year" -> ToExpression[StringTake[#Code, -4]]|> &];

is equivalent to

tdata = tdata[
  All, <|#, "Region" -> ToExpression[StringTake[#Code, 3]], 
    "Year" -> ToExpression[StringTake[#Code, -4]]|> &]

Another way to do the mapping from Region number to State entity

stateRegionMap = EntityList[
  EntityClass["AdministrativeDivision", "ContinentalUSStates"]] // 
 AssociationThread[Range[Length[#]], #] &    

tdata = tdata[All, <|#, "StateName" -> stateRegionMap[#Region]|> &]
POSTED BY: Rohit Namjoshi
J. M.
J. M.
Posted 6 years ago

I would be very interested to know if anyone can suggest faster, more compact or more readable methods of doing what I am doing here.

I suppose the following code would qualify for "compact" (tho it can certainly be compacted further), but tho I find it readable, other people might have different opinions:

(* U.S. states as Entity[] objects *)
continentalStates = 
           EntityList[EntityClass["AdministrativeDivision", "ContinentalUSStates"]];

(* abbreviated month names *)
months = Table[DateString[d, "MonthNameShort"], {d, DateRange[{0, 1}, {0, 12}, "Month"]}];

(* import data as a Dataset and add headers *)
ds = SemanticImport["ftp://ftp.ncdc.noaa.gov/pub/data/cirs/climdiv/climdiv-tmpcst-v1.0.0-20190408",
                    PadRight[{"String"}, 13, "Real"]];
ds = ds[All, AssociationThread[Prepend[months, "CODE"], Range[13]]];

(* per OP, select data corresponding to the continental states... *)
continentalStateData =
           ds[Select[Between[FromDigits[StringTake[#CODE, 3]], {1, 48}] &], All];

(* ...and then select data up to the past year *)
tempHistorical =
    continentalStateData[Select[FromDigits[StringTake[#CODE, -4]] <
                                DateValue[Now, "Year"] &]]

(* compute the weighted mean of the monthly temperatures;
   note the syntax for LeapYearQ[]! *)
meanTempHistorical = tempHistorical[All, <|"State" :> (FromDigits[StringTake[#CODE, 3]] &), 
                                           "MeanAnnualTemperature" :> 
         With[{ml = Slot /@ months}, 
         Function[Mean[WeightedData[ml,
                  {31, 28 + Boole[LeapYearQ[{FromDigits[StringTake[#CODE, -4]]}]], 
                   31, 30, 31, 30, 31, 31, 30, 31, 30, 31}]]]]|>];

(* moving average of mean annual temperatures *)
gwData = KeyMap[Indexed[continentalStates, #] &,
                Normal[meanTempHistorical[GroupBy["State"],
                With[{ma = MovingAverage[#, 10]}, Indexed[ma, -1] - Indexed[ma, -100]] &,
                     "MeanAnnualTemperature"]]];

GeoRegionValuePlot[gwData, PlotLabel -> "Temperature Rise Over 100 Years", 
                   PlotLegends -> SwatchLegend[Automatic,
                      {"0 °C", "1 °C", "2 °C", "3 °C"}, LegendFunction -> "Frame",
                      LegendLabel -> "Temp, °C"]]

I'll let someone else post the resulting picture, since my version of gedanken Mathematica does not have picture-rendering ability.
POSTED BY: J. M.

The figure above is lacking the legend, which actually is a separate graphic: legend
POSTED BY: John Shonder

Thanks Tim, this helps a lot!
POSTED BY: John Shonder
Anonymous User
Anonymous User
Posted 6 years ago

I would be very interested to know if anyone can suggest faster, more compact or more readable methods of doing what I am doing here.

Not if you value your time in making such an alternate thing.

Maybe you should specify Hz if you need a framerate speed you are trying to reach.

You may experience delays with downloading map data from the wolfram server ... the answer there is to look into getting your own copy or speeding that up. It's pretty new how that all works and wolfram has been investigating ways to get the best of both worlds of convenience and options.
POSTED BY: Anonymous User

Thank you for those improvements Rohit.
POSTED BY: John Shonder

Thanks to all who have contributed to improving this. Here's what the map looks like:

US Climate Change
POSTED BY: John Shonder

This representation implies that the climate obeys state boundaries, which is nonsense. I think a truer picture would result from creating a list of the centroids of the states, associating the temperature change data to these centroids, then create an Interpolation function bounded by the outer boundary of the US, then make a DensityPlot map of this interpolation, and then layer on the state boundaries if you must. You would hopefully eliminate goofy things like the jagged temperature change boundary across the central states.
POSTED BY: Christopher Purcell

A public display of Eating your own dog food can be educational. Starting from the post by J.M. above, we can extract the Legend colours from the FullForm of the final plot. We test these as below:

SwatchLegend[{RGBColor[1.`, 0.92`, 0.5`], 
  RGBColor[1.`, 0.6688982203952127`, 0.24889822039521264`], 
  RGBColor[0.8933132970930611`, 0.34531930796514276`, 0.`] , 
  RGBColor[0.7`, 0.21`, 0.`]}, {"0 \[Degree]C", "1 \[Degree]C", 
  "2 \[Degree]C", "3 \[Degree]C"}, LegendFunction -> "Frame", 
 LegendLabel -> "Temp, \[Degree]C"]

We will need the bounding box for the continental states.

{{latmin, latmax}, {lonmin, lonmax}} = 
 GeoBounds[
  EntityClass["AdministrativeDivision", "ContinentalUSStates"]]

Next add the state centroids to the Association. The list continentalStates was defined in the post by J.M. above.

Map[(gwData[#] = 
    Flatten[{ GeoPosition[#][[1]], 
      gwData[#]}]) &, continentalStates] ;    (* remap the gwData 
Association so each element contains the centroid's latitude,longitude,temperature. *)

 GeoGraphics[ GeoMarker /@ Normal[gwData ][[All, 2, {1, 2}]]  ,  
  GeoRange -> {{latmin, latmax}, {lonmin, lonmax}}, 
 GeoRangePadding -> None]  (* plot the state centroids *)

{Min[Normal[gwData ][[All, 2]]  [[All, 3]]]
 , Max[Normal[gwData ][[All, 2]]  [[All, 
    3]]]} (* this is the range of temperature changes in the data *)

We have temps outside the range used in the above plot's Legend. So we may have to extend our Legend because DensityPlot may clip. We add another legend value.

SwatchLegend[{RGBColor[1.`, 0.92`, 0.5`], 
  RGBColor[1.`, 0.6688982203952127`, 0.24889822039521264`], 
  RGBColor[0.8933132970930611`, 0.34531930796514276`, 0.`] , 
  RGBColor[0.7`, 0.21`, 0.`], 
  RGBColor[0.5`, 0.1`, 0.`]}, {"0 \[Degree]C", "1 \[Degree]C", 
  "2 \[Degree]C", "3 \[Degree]C", "4 \[Degree]C"}, 
 LegendFunction -> "Frame", LegendLabel -> "Temp, \[Degree]C"]

We create an Interpolation function f that accepts {long, lat} values and returns temperatures.

f = Interpolation[
  Map[{{#[[2]], #[[1]]}, #[[3]]} &, Normal[gwData ][[All, 2]] ], 
  InterpolationOrder -> 
   1] (* we have to reverse {lat,long} order in argument list to get \
sensible alignment of the plot *)

region =   
 Rectangle[{lonmin, latmin}, {lonmax, 
   latmax} ] ; (* define a Region that includes all the continental \
states. Use this to define the DensityPlot *)

   plot = DensityPlot[ f[u, v] , {u, v} \[Element] region, 
     ColorFunctionScaling -> False, 
     ColorFunction -> (Blend[ {{0, RGBColor[1.`, 0.92`, 0.5`]}, {1, 
           RGBColor[1.`, 0.6688982203952127`, 0.24889822039521264`]}, {2,
            RGBColor[0.8933132970930611`, 0.34531930796514276`, 
            0.`]} , {3, RGBColor[0.7`, 0.21`, 0.`]},
          {4, RGBColor[0.5`, 0.1`, 0.`]}}, #] &),  
      Frame -> False, AspectRatio -> 1, PlotRangePadding -> None, 
     ImagePadding -> None, PlotRange -> All]  (* remove labels and axes *)

Rasterize this plot and insert it as a GeoImage into a GeoGraphics.

 pp = GeoGraphics[{GeoStyling[{"GeoImage",  
      Rasterize[  plot, RasterSize -> 200]}], EdgeForm[Thin], 
    Polygon[Entity["Country", "UnitedStates"]]}, 
   PlotLabel -> "Temperature Rise Over 100 Years",
   Epilog -> {Inset[
      SwatchLegend[{RGBColor[1.`, 0.92`, 0.5`], 
        RGBColor[1.`, 0.6688982203952127`, 0.24889822039521264`], 
        RGBColor[0.8933132970930611`, 0.34531930796514276`, 0.`] , 
        RGBColor[0.7`, 0.21`, 0.`],
        RGBColor[0.5`, 0.1`, 0.`]}, {"0 \[Degree]C", "1 \[Degree]C", 
        "2 \[Degree]C", "3 \[Degree]C", "4 \[Degree]C"}, 
       LegendFunction -> "Frame", 
       LegendLabel -> "Temp, \[Degree]C"] , {Right, Bottom}, {Right, 
       Bottom}]}]  ;
states = GeoGraphics[ {GeoStyling[Opacity[0]], EdgeForm[Thin], 
    Polygon[EntityClass["AdministrativeDivision", 
      "ContinentalUSStates"]]}, GeoBackground -> "CountryBorders"]  ;
Show[pp, states] (* add state boundaries to the Density plot *)

Interpolated Temperature Rise

I do not fully understand GeoStyling[{"GeoImage", ...}] so am guessing that the Interpolation function has invented the hot spot on the Texas Gulf coast because extrapolation was used (south and east of the centroid). In California, I don' t know what happened, except California appears cooler than it should be.

So although very pretty, this map is even more nonsensical than the original (where climate obeyed state boundaries), and demonstrates that Interpolation outside the convex hull of the state centroids is dangerous. We need more data points.
POSTED BY: Christopher Purcell

I don't know if anyone is trying to use the geospatial statistics method? For example, IDW(Inverse Distance Weighted) or Kriging. I use an example obtained on the Internet to try to make the IDW method.

Attachments:
IDW.nb
POSTED BY: Tsai Ming-Chou
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract