# US Climate Change at the County Level

Posted 2 years ago
10830 Views
|
21 Replies
|
30 Total Likes
|
 (Edited 4/29/2019) Continuing my efforts to visualize climate change in the US, I found that the US National Oceanic and Atmospheric Administration (NOAA) maintains a dataset of monthly average temperatures for every county in the 48 contiguous US States from 1895 to the present. To avoid burdening the server during testing, I downloaded the file to my computer, but it is available here. On the NOAA ftp site the file will have the name climdiv-tmpccy-v1.0.0-YYYYMMDD, where YYYY, MM and DD correspond to the year, month and date when the file was updated. When I downloaded the file its name was climdiv-tmpccy-v1.0.0-20190408. It takes some time to read in all 388,375 lines of this file: cdata = SemanticImport["climdiv-tmpccy-v1.0.0-20190408.txt", PadRight[{"String"}, 13, "Real"]]; The dataset consists of 13 columns: the first is a string of digits encoding the year and something that is almost -- but not quite -- the FIPS code for the county. The next 12 are the monthly average temperatures in degrees Fahrenheit in that county for January through December of that year. For clarity, I like to add column headers first: cdata = cdata[All, AssociationThread[ Prepend[Table[ DateString[d, "MonthNameShort"], {d, DateRange[{0, 1}, {0, 12}, "Month"]}], "Code"], Range[13]]]; To correct the FIPS codes, we require the FIPS codes for the 48 contiguous States: stateFIPS = EntityValue[ EntityList[ EntityClass["AdministrativeDivision", "ContinentalUSStates"]], "FIPSCode"]; Now we can add a column containing the correct FIPS code for each county. While we're at it we add a column for the year: cdata = cdata[ All, <|"FIPS" -> stateFIPS[[ToExpression[StringTake[#Code, 2]]]] <> StringTake[#Code, {3, 5}], "Year" -> ToExpression[StringTake[#Code, -4]], #|> &]; Then I build an association of FIPS codes and entity values for the corresponding counties: counties = EntityList[ EntityClass["AdministrativeDivision", "USCountiesAllCounties"]]; massoc = AssociationThread[EntityValue[counties, "FIPSCode"], counties]; With this association, we can add a column to the dataset corresponding to the entity value of each county. In the same statement I also add a column for the average annual temperature: cdata = cdata[ All, <|"County" -> massoc[#FIPS], #, "Tavg" -> Mean[WeightedData[{#Jan, #Feb, #Mar, #Apr, #May, #Jun, #Jul, \ #Aug, #Sep, #Oct, #Nov, #Dec}, {31, 28 + Boole[LeapYearQ[#Year]], 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}]]|> &]; The file includes data for the first couple of months of 2019, so I strip that out: cdata = cdata[Select[#Year <= 2018 &]]; I tried several different ways of smoothing the temperatures and decided on an moving average of the past 10 years: s = cdata[GroupBy["County"], With[{ma = MovingAverage[#, 10]}, Last[ma] - ma[[-100]]] &, "Tavg"]; I experimented with the color functions available in Mathematica, but decided to build my own using a color pallette from the ColorBrewer website. col1 = RGBColor[{255, 245, 240}/255]; col2 = RGBColor[{254, 224, 210}/255]; col3 = RGBColor[{252, 187, 161}/255]; col4 = RGBColor[{252, 146, 114}/255]; col5 = RGBColor[{251, 106, 74}/255]; col6 = RGBColor[{239, 59, 44}/255]; col7 = RGBColor[{203, 24, 29}/255]; col8 = RGBColor[{165, 15, 21}/255]; col9 = RGBColor[{103, 0, 13}/255]; qq = Values[s]; zz = Table[Min[qq] + i*(Max[qq] - Min[qq])/9, {i, 9}]; cfunc[x_?NumericQ] := Which[ x <= zz[[1]], col1, x <= zz[[2]], col2, x <= zz[[3]], col3, x <= zz[[4]], col4, x <= zz[[5]], col5, x <= zz[[6]], col6, x <= zz[[7]], col7, x <= zz[[8]], col8, x <= zz[[9]], col9] Then it's just a matter of building the legend and displaying the map: legend = SwatchLegend[{col1, col2, col3, col4, col5, col6, col7, col8, col9}, {"-0.30 - 0.28", " 0.28 - 0.87", " 0.87 - 1.45", " 1.45 - 2.04", " 2.04 - 2.62", " 2.62 - 3.21", " 3.21 - 3.79", " 3.79 - 4.38", " 4.38 - 4.96"}, LegendMarkers -> Graphics[{EdgeForm[Black], Opacity[1], Rectangle[]}], LegendLabel -> "\[CapitalDelta]T(\[Degree]F)", LegendFunction -> (Framed[#, RoundingRadius -> 5] &), LegendMargins -> 5]; climvis = GeoRegionValuePlot[s, Frame -> True, FrameTicks -> None, FrameLabel -> {"Change in Mean Annual Temperature for US Counties, \ 1919-2018"}, LabelStyle -> Larger, PlotLegends -> Placed[legend, Right], ColorFunction -> cfunc, ColorFunctionScaling -> False, PlotStyle -> Directive[EdgeForm[{Thin, White}]], GeoBackground -> None, GeoProjection -> {"LambertAzimuthal", "Centering" -> GeoPosition[{30, -195/2}]}, PlotRange -> {{-0.37, 0.38}, {-0.13, 0.38}}, ImageSize -> 800, PlotLegends -> Placed[legend, Right]] And finally here is the result. This is of course just a "proof of concept" exercise. There are any number of ways to analyze this dataset and visualize the results. Attachments:
21 Replies
Sort By:
Posted 2 years ago
 Hi John, Thanks for sharing. It seems a great code. But unfortunately I don't have access to the website and data (see screenshot below). Please can you support me with the data? Your support will be highly appreciated ! Best Regards,....Jos
Posted 2 years ago
Posted 2 years ago
 Hi John, and thanks for sharing - very educational. One thing I don't fully understand (and maybe you can explain it a bit?) is how you calculated the values for the legend.E.g. when I'm retracing your steps, I don't see any negative values in qq: Min[qq] 0.134336 and zz {0.74857, 1.3628, 1.97704, 2.59127, 3.20551, 3.81974, 4.43397, \ 5.04821, 5.66244} What am I missing?
Posted 2 years ago
 My mistake. I originally used a moving average of 10 years for the smoothing, and then changed it to an exponential average. The values in the legend correspond to the moving average case. Looking at both figures, there really is not much of a difference. This is a problem when using a bespoke legend as I did. I'm sure it would be possible to generate the legend automatically with the correct limits -- and perhaps someone will be interested enough to show us how to do that -- but for now I went back and edited my original post to correct everything. The figure and the file now use the moving average.
Posted 2 years ago
 - Congratulations! This post is now featured in our Staff Pick column as distinguished by a badge on your profile of a Featured Contributor! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!
Posted 2 years ago
 Excellent analysis! A quick question, what was the standard deviation in temperature measurements? I am asking that because some Delta T are small and I wonder to what extend we can interpret them as a meaningful shift in T?
Posted 2 years ago
 Mads, that is a great question. All of the data I used are online at the site I referenced, and it would be easy to modify my code to get a map of the standard deviation by county.
Posted 2 years ago
 John, I may be mistaken in decoding the file, but I believe the file shows about 17 county entries for Michigan and it has 83. For Arizona, the state code is 04 and has 15 counties. The file shows 58 counties being reported for 04. FYI. Thanks for the information as well as your analysis. Very Interesting.
Posted 2 years ago
 Ed, you have to read the documentation for that file. As I said above, it does not contain the FIPS code, but a number that can be transformed into a FIPS code. And I showed how to do it.
Posted 2 years ago
 Thank you. Must be my error in calculating the FIPS.-Ed
Posted 2 years ago
 John, the documentation's state code reference was the key for me. Thank you again for the suggestion.-Ed
Posted 2 years ago
 The question I thought someone was going to ask is where did they get temperature data for counties that didnt even exist at the time. Arizona for example didnt become a state until 1912 and yet there are the average monthly temperatures for each of its current counties back to 1895. The answer is through kriging.
Posted 2 years ago
 Great chart with potential for display of any year span.
Posted 2 years ago
 John - I used your post for a project with my kids this morning. Thanks for sharing this amazing work - they loved it and had a few ideas that they'd like to explore a bit more. Can't wait to work on these ideas with them. Here's what we did:https://mikesmathpage.wordpress.com/2019/06/16/using-john-shonders-amazing-us-temperature-visualization-with-kids/
Posted 2 years ago
 Great way to get the kids interested. I did much the same as you but experimented with different start years. TBP.BTW the missing county is entirely an Indian reservation with no county seat. (I wonder why kriging was not used).
Posted 2 years ago
 Mike, how nice to know you were able to use the code and data to teach your children. That's why I posted it: hoping people would extend the analysis, try out different smoothing techniques, use different starting and ending points, etc. I watched your videos and it sounds like your boys have some good ideas. I'll be interested to see what they come up with. Happy Father's Day.
Posted 2 years ago
 I decided to publish my changes. I separated the temperature ranges at 0 so there's no doubt as to the direction of change. Also the start date is 20 years later to see how the trend changes.
Posted 2 years ago
 Very interesting Doug! In later versions of this map I too decided to use shades of blue for negative changes. It makes sense visually. You also see that the starting point makes a difference too. The map at the top of this page shows temperature changes between 1991 and 2012 compared to the average temperature between 1901 and 1960.
Posted 2 years ago
 Starting point makes a difference - 1991 is a unique starting point. Mt. Pinatubo eruption. Nearly 20 million tons of sulfur dioxide were injected into the stratosphere in Pinatubo's 1991 eruptions, and dispersal of this gas cloud around the world caused global temperatures to drop temporarily (1991 through 1993) by about 1°F (0.5°C). (https://pubs.usgs.gov/fs/1997/fs113-97/) I'll attach a run of your 10-year moving averages 1991-2012, not as year-specific but interesting. (Is there a way to embed the legend in the chart?)!legend for climate by counties][1] Attachments: