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

Posted 1 year ago
10465 Views
|
12 Replies
|
19 Total Likes
|
 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}}] 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}]] 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]]; 
12 Replies
Sort By:
Posted 1 year ago
 -- you have earned Featured Contributor Badge 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 1 year ago
 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}]&]. 
Posted 1 year ago
 Alex, I tried the code in both V12.1 and V12.0 released version and the piece of code works.
Posted 1 year ago
 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 1 year ago
 Does it work if you change the first line to this? COVID19CountyData = ResourceFunction["NYTimesCOVID19Data"]["USCountiesTimeSeries"]; 
Posted 1 year ago
 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. 
Posted 1 year ago
 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.
Posted 1 year ago
 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}}] 
Posted 1 year ago
 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"] 
Posted 1 year 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.
 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?