Who has not encountered a stink bug? Perhaps the better question is not if, but when. I remember well my first interactions with stink bugs—partly because of their pungent, cilantro-like odor, but also because in my native Catalan language they are called Bernat pudent (“stinky Bernat”) and Bernat is my twin brother’s name.

So when I encountered the stink bug again when visiting Champaign, Illinois, for the 2019 Wolfram Technology Conference, it brought up a lot of fond childhood memories. This time, however, two things had changed: the frequency of encounters with the stink bug seemed exponentially greater, and I now had the Wolfram Language to more fully (and computationally) satisfy my curiosity about this reviled insect and its growing impact on our ecosystem. So to get a better picture of the arrival and spread of this invasive bug across the US, I used available observation data and the Wolfram Language to make a map of sightings over the past two decades.

Importing Distribution Data

In the US, there are dozens of native stink bug species, but in this post I’m going to focus mainly on a non-native species that has become quite common in North America and fairly unpopular over the last two decades: the brown marmorated stink bug, also known by its scientific name as Halyomorpha halys (Stål, 1855). To get as comprehensive a dataset as possible, we will be importing stink bug data from three sources. After saving the data files in the same folder as our notebook, we should set the directory path appropriately:

Engage with the code in this post by downloading the Wolfram Notebook

&#10005 SetDirectory[NotebookDirectory[]];

We’ll start by importing data on stink bug distribution from the EPPO Global Database. Using SemanticImport, we can specify data interpretations for the desired columns and get a Dataset object for easy exploration. In this case, Interpreter["Country"] is being passed as a pure function for richer interpretation:

&#10005 dataEPPO = SemanticImport[ "distribution_HALYHA.csv", <| "country code" -> Interpreter["Country"], "state" -> "String", "Status" -> "String"|>, "HeaderLines" -> 1]

We can get a quick overview by using Select to get national data from around the world (represented by items with a blank “state” property):

&#10005 nationalData = dataEPPO[Select[StringMatchQ[#state, ""] &]];

The data assigns each sighting to one of six distinct status levels—either absent or present, and the quality/intensity of the report:

&#10005 Column[Union@Normal@nationalData[[All, "Status"]]]

We can use GroupBy to collect the countries by status:

&#10005 countriesByStatus = KeySort@nationalData[GroupBy[#Status &]][[All, All, "country code"]] // Normal;

Using GeoRegionValuePlot, we can then create a map of stink bug distribution and level of occurrence per country, based on the EPPO data:

&#10005 GeoRegionValuePlot[ Thread /@ Thread[Values@countriesByStatus -> {RGBColor[ 0.00784313725490196, 0.5098039215686274, 0.9294117647058824], RGBColor[ 0.5411764705882353, 0.7137254901960784, 0.027450980392156862`], RGBColor[ 0.9333333333333333, 0.5098039215686274, 0.9333333333333333], RGBColor[ 0.996078431372549, 0.9882352941176471, 0.03529411764705882], RGBColor[ 0.996078431372549, 0.3607843137254902, 0.027450980392156862`], RGBColor[1, 0, 0]}] // Flatten, { PlotStyle -> Opacity[0.6], ImageSize -> 550, GeoBackground -> "CountryBorders", Epilog -> Inset[ Framed[ SwatchLegend[{ RGBColor[0.00784313725490196, 0.5098039215686274, 0.9294117647058824], RGBColor[0.5411764705882353, 0.7137254901960784, 0.027450980392156862`], RGBColor[0.9333333333333333, 0.5098039215686274, 0.9333333333333333], RGBColor[0.996078431372549, 0.9882352941176471, 0.03529411764705882], RGBColor[0.996078431372549, 0.3607843137254902, 0.027450980392156862`], RGBColor[1, 0, 0]}, { "Absent, intercepted only", "Absent, unreliable record", "Present, few occurrences", "Present, no details", "Present, restricted distribution", "Present, widespread"}, LegendLayout -> "Row"], RoundingRadius -> 5, Background -> Directive[ GrayLevel[1], Opacity[0.9]]], {35, -63}]}]

Now let’s import the data from iNaturalist and EDDMapS (Early Detection and Distribution Mapping System), both of which contain crowdsourced observation data. We’ll use SemanticImport again. Because the locations are listed as latitude/longitude pairs, no custom interpretation is needed:

&#10005 dataEDDMapS = SemanticImport["EDDMapS_Halys.csv"]; dataiNaturalist = SemanticImport["iNaturalist_US_Stink_bugs.csv"];

The iNaturalist data contains observations of many stink bug species; for this exploration, we only need observations of the invasive stink bug, Halyomorpha halys:

&#10005 dataiNaturalist = dataiNaturalist@ Select[#["scientific_name"] == "Halyomorpha halys" &];

With our data imported, we’re now ready to create a distribution map.

Mapping Out the Invasion

The main purpose of this visualization is to determine geographic distribution throughout time. I decided the easiest representation for this would be an interactive, color-coded map showing the combined data across several years. First, we assign a color to each year since 1998, roughly when stink bugs were first sighted:

&#10005 colorYear[year_Integer] := Part[Table[ ColorData[{"SolarColors", "Reverse"}][i], {i, 0, 1, 0.05}], year - 1998];

Then we need to extract coordinates. Since the datasets have different formats, we need separate functions for selecting the observations from each. For EDDMapS, we group by observation date, extracting the latitude and longitude from each data point:

&#10005 EDDMapSCoords = Values /@ Normal[dataEDDMapS@GroupBy[#["ObsDate"]["Year"] &]][[All, All, {"Latitude", "Longitude"}]];

The process is the same for the iNaturalist dataset, but with different key names:

&#10005 iNaturalistCoordinates = Values /@ Normal[dataiNaturalist@GroupBy[#["observed_on"]["Year"] &]][[All, All, {"latitude", "longitude"}]];

We can use Merge and DeleteMissing to combine and clean the lists:

&#10005 combinedCoords = Replace[DeleteMissing[Flatten[#, 1], 2], {} -> Nothing, 3] & /@ Merge[{EDDMapSCoordinates, iNaturalistCoordinates}, Identity];

Then we create an array of geodetic positions with GeoPosition:

&#10005 positionsYear = GeoPosition /@ combinedCoords;

Here’s a styled background for our animation, using GeoListPlot to display polygons for each state over a satellite image:

&#10005 bgImage = GeoListPlot[ EntityList@ EntityClass["AdministrativeDivision", "AllUSStatesPlusDC"], { PlotStyle -> Directive[ EdgeForm[{ GrayLevel[0.85], Thickness[Tiny], Opacity[0.8]}], FaceForm[{ Opacity[0.2], GrayLevel[1]}]], GeoBackground -> GeoStyling["Satellite", Opacity[0.8]], GeoRange -> Entity["Country", "UnitedStates"], PlotLegends -> None, ImageSize -> 550}]

We also want to include a thumbnail image of the stink bug on the final GIF. The following specimen was trying to overwinter inside my flat in Barcelona (we also have the same invasive stink bug species in Europe):

&#10005 thumbnail = CloudGet["https://wolfr.am/JNu0LzR8"];

Finally, we can combine my thumbnail image, the custom color function and the previous positions-per-year function in GeoListPlot (with other styling tweaks enclosed in iconized code) to display successive sightings per year starting in 2001, using Show to superimpose the data onto our background:

&#10005 distributionYear[year_Integer] := Show[bgImage, GeoListPlot[ Table[positionsYear[y], {y, 2001, year}], ({ PlotMarkers -> Point, PlotStyle -> Table[ Directive[ colorYear[y], Opacity[0.5], PointSize[0.0045]], {y, 2001, #}], GeoRange -> Entity["Country", "UnitedStates"], PlotLegends -> None, ImageSize -> 550, Epilog -> { Inset[ Framed[ Style[ ToString[#], colorYear[#], FontSize -> 36, FontFamily -> "Helvetica"], Background -> Directive[White, Opacity[0.2]], FrameStyle -> Transparent, RoundingRadius -> 5], Scaled[{0.5, 0.5}]], Inset[ Column[{ ImageResize[thumbnail, 75], Style["Halyomorpha halys", Italic, Bold, FontSize -> 13, LightGray]}, Alignment -> Right], Scaled[{0.86, 0.23}]]}}& )[year]], GeoRange -> Entity["Country", "UnitedStates"]]

The resulting map shows the accumulation of data points up to a given year, along with a central label showing the year:

&#10005 distributionYear[2017]

Now that we can generate distribution maps of each year, it is quite easy to assemble them into an animated GIF. We’re essentially creating a ParallelTable of the different years’ maps, then exporting as an animation (repeating the final frame for a longer viewing time):

&#10005 Export["Halys_US_Distribution.gif", ParallelTable[ distributionYear[y], {y, Append[Range[2001, 2019], 2019]}], "DisplayDurations" -> 0.6, "AnimationRepetitions" -> Infinity]

Seeing the result assembled into an animation, it’s immediately apparent that the reported stink bug sightings have increased dramatically and over an increasing amount of territory as the years continue. Over 20 years, there is seemingly unrestricted growth—likely due to the absence of natural predators like the samurai wasp, which help control stink bug populations in their native habitat in East Asia. And thanks to the contributions of citizen scientists on platforms like iNaturalist, we are able to use their data to map these sightings more accurately.

When I posted my GIF to Reddit it touched a nerve, with more than 41,000 Reddit users engaging with the post! And circling back to my twin brother Bernat: he had shown the GIF to some friends in Providence, Rhode Island, where he was living back then—and two of them told him that they had already seen the GIF on the front page of Reddit (demonstrating a different kind of “invasive” distribution—internet virality!).

The Weight of Imbalance

The distribution map has provided an easily comprehensible confirmation of a shared experience: the invasive brown marmorated stink bug has spent the past two decades becoming a near-ubiquitous nuisance. But even more significantly, we can see that they are a serious threat in the “breadbasket” states (like Illinois), where a majority of the country’s food crops are grown. The West Coast, geographically opposite the location of the stink bug’s original arrival, has also developed its own share of problems with the insect. Even California—where orchards and crops are already threatened by fire and drought—isn’t safe. Accounts of an independent West Coast arrival event (confirmed by genetic testing) are clearly borne out by the distribution map data.

We may be irritated by the one, five or even some thousands of stink bugs we encounter ourselves, but seeing 20 years’ worth of an invasive species’ unabated spread across the nation better speaks to the full weight of imbalanced biodiversity. I encourage you to sign up for crowdsourced sites like iNaturalist to become a citizen scientist, and start contributing to science with your own observations. As more people engage in the conversation, the better chance we have to stop a 30-year retrospective on the stink bug invasion!

Join us—share your own computational explorations or stink bug observations via the Wolfram Community thread that I initially started, or check out my other posts on biodiversity. Every data point is important!

Additional Resources

“Brown Marmorated Stink Bug.” National Invasive Species Information Center, USDA.
  www.invasivespeciesinfo.gov/profile/brown-marmorated-stink-bug.

Rice, K. B., C. J. Bergh, E. J. Bergmann, et al. “Biology, Ecology, and Management of
  Brown Marmorated Stink Bug (Hemiptera: Pentatomidae).” Journal of Integrated
  Pest Management 5.3 (2014): A1–A13.