Message Boards Message Boards

[GIF] Distance to nearest confirmed US COVID-19 case

Posted 5 years ago

enter image description here

The New York Times recently published COVID-19 time series data for each US county. This gives us a very granular look at the spread of the illness over time. This will let us approximate the distance to the nearest confirmed COVID-19 case for any given location in the US.


To start this task, first we import the data, canonicalize each time series to start on Jan 21, and only keep the lower 48 states:

COVID19CountyData = ResourceFunction["NYTimesCOVID19Data"]["USCountiesTimeSeries"];

COVID19CountyData = COVID19CountyData[All, "Cases", 
   TimeSeries[TimeSeriesInsert[#, {DateObject[{2020, 1, 20}], 0}], ResamplingMethod -> {"Interpolation", InterpolationOrder -> 0}] &];

COVID19CountyData = COVID19CountyData[KeySelect[!MissingQ[#] && !MatchQ[#, Entity[_, {_, "Alaska" | "Hawaii", _}]] &]];

Next we'll sample points in the US by discretizing its polygon, and then interpolating distances from these samples:

usstates = EntityClass["AdministrativeDivision", "ContinentalUSStates"];
usa = RegionUnion[BoundaryDiscretizeGraphics /@ usstates["Polygon"][[All, All, 1]]];
usatri = DiscretizeRegion[usa, MaxCellMeasure -> .05];
usasamp = GeoPosition /@ MeshCoordinates[usatri];

Here are the points we'll be sampling from:

MeshRegion[TransformedRegion[usatri, ReflectionTransform[{1, 0}]@*RotationTransform[π/2]], MeshCellStyle -> {0 -> {PointSize[Small], Red}, 1 -> {Thin, Black}}]

enter image description here

Here's the GraphicsComplex we'll be adding VertexColors to:

usagcomp = GraphicsComplex[
  GeoPosition[MeshCoordinates[usatri]],
  MeshCells[usatri, 2, "Multicells" -> True]
];

For each point, to find the GeoDistance to the nearest affected county, we'll approximate by finding the nearest GeoPosition of the county polygons. This will speed up calculations.

countypolys = EntityValue[Normal[Keys[COVID19CountyData]], "Polygon", "EntityAssociation"];

Finally, we have a function that takes in a date, finds all counties that have confirmed cases then, uses their polygons to find distances to each point we chose in the US, then creates a GeoDistance heat map:

COVID19DistanceMap[date_DateObject?DateObjectQ] :=
  Block[{counties, polys, coords, dists, gcomp, map, legend},
    counties = COVID19CountyData[Select[#[date] > 0&]];

    polys = Lookup[countypolys, Normal[Keys[counties]]];
    coords = GeoPosition /@ Join @@ Cases[polys, _List?MatrixQ, ∞];

    coords = PositionIndex[Counts[coords]][1];
    dists = Clip[QuantityMagnitude[Nearest[coords -> "Distance", usasamp][[All, 1]], "Miles"], {0., 100.}];

    gcomp = Append[usagcomp, VertexColors -> ColorData[{"Rainbow", {-100, 0}}] /@ Minus[dists]];

    map = GeoGraphics[
      {
        GeoStyling[Opacity[1]], gcomp, 
        {Gray, polys},
        {EdgeForm[Black], FaceForm[], Polygon[usstates]}
      },
      ImageSize -> 1024
    ];

    legend = makeLegend[date];

    Legended[map, Placed[legend, {0.1725, 0.125}]]
]

makeLegend[date_] :=
  Framed[
    Column[{
      BarLegend[{ColorData[{"Rainbow", "Reversed"}][.01#]&, {0, 100}},
        LegendLayout -> "Row", ImageSize -> 300, LegendLabel -> Style["\[ThinSpace]   Nearest Confirmed COVID\[Hyphen]19 Case (mi), " <> DateString[date, {"Month","/","Day","/","YearShort"}], 13, GrayLevel[0.2]]],
      Row[{Spacer[24], 
        SwatchLegend[{Gray, ColorData["Rainbow", 0]}, {"Affected area", "> 100"}, LegendMarkerSize -> 13, LabelStyle -> GrayLevel[0.2], LegendLayout -> "Row", Spacings -> {0.5, 10.5}]}]
    },
    Spacings -> 0],
  Background -> GrayLevel[0.95],
  RoundingRadius -> 10,
  FrameStyle -> {AbsoluteThickness[1.5], GrayLevel[0.75]},
  FrameMargins -> {{0, 10}, {10, 15}}
]

Here's the result for March 16, 2020. Notice that we clip the distance at 100 miles. This will give a consistent scale over time.

COVID19DistanceMap[DateObject[{2020, 3, 16}]]

enter image description here

We can now create a graphic for each day and export as a gif. Note that this is somewhat time consuming, but multiple kernels help.

frames = ParallelMap[
  Rasterize[COVID19DistanceMap[#]] &, 
  DateRange[DateObject[{2020, 1, 21}], DateObject[{2020, 3, 24}]], 
  Method -> "FinestGrained"
];

Export["covid19_distance_map.gif", frames, 
  AnimationRepetitions -> ∞, 
  "DisplayDurations" -> Append[ConstantArray[.2, Length[frames] - 1], 3]];
POSTED BY: Greg Hurst
12 Replies
Posted 4 years ago

Wow, this is really amazing to see this in a gif format. Great job! Thanks for sharing this!

POSTED BY: Yasmin Hussain

I'm on 12.1, and my kernel crashes whenever I try to execute

usatri = DiscretizeRegion[usa, MaxCellMeasure -> .05];

Anyone else having this problem? Maybe it's lack of RAM?

POSTED BY: Andrew Dabrowski

Thank you. Finally, I reproduced the result by changing the code as follows

COVID19CountyData = 
  ResourceFunction["NYTimesCOVID19Data"]["USCountiesTimeSeries"];

COVID19CountyData = 
  COVID19CountyData[All, "Cases", 
   TimeSeries[TimeSeriesInsert[#, {DateObject[{2020, 1, 20}], 0}], 
     ResamplingMethod -> {"Interpolation", 
       InterpolationOrder -> 0}] &];

COVID19CountyData = 
  COVID19CountyData[
   KeySelect[! MissingQ[#] && ! 
       MatchQ[#, Entity[_, {_, "Alaska" | "Hawaii", _}]] &]];

coord = CountryData["USA", "Polygon"][[1, 1, 1]];
bmr = BoundaryMeshRegion[coord, Line[Append[Range[Length[coord]], 1]]];

usatri = TriangulateMesh[bmr, MaxCellMeasure -> .1];

MeshRegion[
 TransformedRegion[usatri, 
  ReflectionTransform[{1, 0}]@*RotationTransform[\[Pi]/2]], 
 MeshCellStyle -> {0 -> {PointSize[Small], Red}, 1 -> {Thin, Black}}]
usagcomp = 
  GraphicsComplex[GeoPosition[MeshCoordinates[usatri]], 
   MeshCells[usatri, 2, "Multicells" -> True]];

usstates = 
  EntityClass["AdministrativeDivision", "ContinentalUSStates"];

usasamp = GeoPosition /@ MeshCoordinates[usatri];


countypolys = 
  EntityValue[Normal[Keys[COVID19CountyData]], "Polygon", 
   "EntityAssociation"];
COVID19DistanceMap[date_DateObject?DateObjectQ] := 
 Block[{counties, polys, coords, dists, gcomp, map, legend}, 
  counties = COVID19CountyData[Select[#[date] > 0 &]];
  polys = Lookup[countypolys, Normal[Keys[counties]]];
  coords = 
   GeoPosition /@ Join @@ Cases[polys, _List?MatrixQ, \[Infinity]];
  dists = 
   Clip[QuantityMagnitude[
     Nearest[coords -> "Distance", usasamp][[All, 1]], "Miles"], {0., 
     100.}];
  gcomp = 
   Append[usagcomp, 
    VertexColors -> ColorData[{"Rainbow", {-100, 0}}] /@ Minus[dists]];
  map = GeoGraphics[{GeoStyling[Opacity[1]], 
     gcomp, {Gray, polys}, {EdgeForm[Black], FaceForm[], 
      Polygon[usstates]}}, ImageSize -> 1024];
  legend = makeLegend[date];
  Legended[map, Placed[legend, {0.1725, 0.125}]]]

makeLegend[date_] := 
 Framed[Column[{BarLegend[{ColorData[{"Rainbow", 
          "Reversed"}][.01 #] &, {0, 100}}, LegendLayout -> "Row", 
     ImageSize -> 300, 
     LegendLabel -> 
      Style["\[ThinSpace]   Nearest Confirmed COVID\[Hyphen]19 Case \
(mi), " <> DateString[date, {"Month", "/", "Day", "/", "YearShort"}], 
       13, GrayLevel[0.2]]], 
    Row[{Spacer[24], 
      SwatchLegend[{Gray, ColorData["Rainbow", 0]}, {"Affected area", 
        "> 100"}, LegendMarkerSize -> 13, 
       LabelStyle -> GrayLevel[0.2], LegendLayout -> "Row", 
       Spacings -> {0.5, 10.5}]}]}, Spacings -> 0], 
  Background -> GrayLevel[0.95], RoundingRadius -> 10, 
  FrameStyle -> {AbsoluteThickness[1.5], GrayLevel[0.75]}, 
  FrameMargins -> {{0, 10}, {10, 15}}]

But the last piece of visualization code is very slow. How can this be accelerated?

frames = ParallelMap[Rasterize[COVID19DistanceMap[#]] &, 
  DateRange[DateObject[{2020, 3, 21}], DateObject[{2020, 3, 25}]], 
  Method -> "FinestGrained"]

Figure 1

Posted 5 years ago

Yes, the slowness comes from the call to Nearest. It's using the GeoDistance metric on 100,000's of points as my way of approximating distance to a polygon with a non-Euclidean metric.

See my edit to COVID19DistanceMap though. I added a line which will eliminate about half of the points beforehand, making it twice as fast.

POSTED BY: Greg Hurst

In which version did you test this code? In version 12.1, this does not work. There is a message

$RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of COVID19CountyData[All,Cases,TimeSeries[TimeSeriesInsert[#1,{DateObject[{2020,1,20}],0}],ResamplingMethod->{Interpolation,InterpolationOrder->0}]&].

Alex, I tried the code in both V12.1 and V12.0 released version and the piece of code works.

work

POSTED BY: Shenghui Yang

Sorry, but it doesn’t work. Perhaps this is a bug of the type discussed on https://mathematica.stackexchange.com/questions/216746/user-interface-mathematica-12-1-terribly-slow

Posted 5 years ago

Does it work if you change the first line to this?

COVID19CountyData = ResourceFunction["NYTimesCOVID19Data"]["USCountiesTimeSeries"];
POSTED BY: Updating Name

Thank you, this piece of code works now. But the second piece of code does not work in v.12.1,

usstates = 
  EntityClass["AdministrativeDivision", "ContinentalUSStates"];
usa = RegionUnion[
   BoundaryDiscretizeGraphics /@ usstates["Polygon"][[All, All, 1]]];
usatri = TriangulateMesh[usa, MaxCellMeasure -> .05];
usasamp = GeoPosition /@ MeshCoordinates[usatri];

We have messages

BoundaryMeshRegion::holm: Unable to find points not in the region for all of the inner boundaries.

BoundaryMeshRegion::holm: Unable to find points not in the region for all of the inner boundaries.

BoundaryMeshRegion::could not make representation: -- Message text not found -- (TriangleLinkMesh)

ToElementMesh::fememins: The mesh elements are not valid. A set of valid mesh element incidents needs to be positive integers and be able to form a complete sequence starting from 1 to the largest incident present. There are missing incidents; a complete sequence cannot be formed.

ToElementMesh::femtemnm: A mesh could not be generated.

I do get the same two BoundaryMeshRegion::holm messages from TriangulateMesh, but I don't get the other messages you see, also using 12.1 (in macOS). The computation succeeds.

I used the standard code to triangulate the mesh

coord = CountryData["USA", "Polygon"][[1, 1, 1]];
bmr = BoundaryMeshRegion[coord, Line[Append[Range[Length[coord]], 1]]];

usatri = TriangulateMesh[bmr, MaxCellMeasure -> .1];

MeshRegion[
 TransformedRegion[usatri, 
  ReflectionTransform[{1, 0}]@*RotationTransform[\[Pi]/2]], 
 MeshCellStyle -> {0 -> {PointSize[Small], Red}, 1 -> {Thin, Black}}]

enter image description here -- you have earned Featured Contributor Badge enter image description here

Your exceptional post has been selected for our editorial column Staff Picks http://wolfr.am/StaffPicks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: EDITORIAL BOARD
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract