Message Boards

GROUPS:

Staff Picks Biological Sciences Data Science Medical Sciences Mathematica Graphics and Visualization Import and Export Wolfram Language Modeling Geographic Information System

Diego Zviovich

Modeling the spatial spread of infection diseases in the US

Diego Zviovich

Posted 2 years ago

9698 Views

5 Replies

14 Total Likes

Follow this post

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus

Introduction

The following document will be the first one on a series to try to replicate the model used by GLEAMviz , reducing the scope to the 48 contiguous US states (US). The focus will be in this region as data for the US is more readily available than in other geographies.
This model

uses socio-demographic and mobility data to simulate the spread of epidemics across a given geography. We’ll explore the spatio-temporal pattern of an epidemic spreading by seeding the starting locus at the different main points of entry.

Model

The model is composed of three different data layers: the population layer, the mobility layer and the disease layer. To build these layers, let’s proceed to capture the available data from the different sources.

The population layer

Our population layer will be built upon the University of Columbia’s “Gridded Population of the World” dataset

. One of the available datasets is the “population count v4.11”

. We’ll download the population data for 2020 at a 30 arc minute resolution (approximately 55km x 55 km). https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-count-rev11/data-download. The population layer will be constituted by a voronoi-like tesellation of the US around the main airports across the country, which we’ll be worked on in the next section. Each subpopulation will be assigned to a cell associated to the closest airport.

Let’s set up a structure to store our work as we move along the project. This structure will be made available in either github or the wolfram data repository. Link will be provided later.

data=AssociationMap[<||>&,{"sources","rawData","temp","results"}]

Out[]=

sources,rawData,temp,results

In[]:=

data["sources","populationData"]=Hyperlink["https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-count-rev11/data-download"];

Data has been saved in the local drive. We’ll proceed to read the data.

In[]:=

data["sources","populationFiles"]=FileNames["*min.asc",FileNameJoin[{NotebookDirectory[],"datasets","population","data"}]];

In[]:=

data["rawData","population"]=Import[First@data["sources","populationFiles"],"Data"];

Let’s review the data structure of the read file. We can observe that it is a matrix of 720 columns and 360 rows, each cell representing a 0.5x0.5 angular degree sector.

data["rawData","population"][[;;6]]

Out[]=

{ncols 720,nrows 360,xllcorner -180,yllcorner -90,cellsize 0.5,NODATA_value -9999}

Processing the data.

In[]:=

data["temp","population"]=Association[];data["temp","population","lines"]=StringReplace[#,{" "",","e"" 10^","-9999""0."}]&/@data["rawData","population"][[7;;]];data["temp","population","numbers"]=ToExpression@StringSplit[#,","]&/@data["temp","population","lines"];

MatrixPlot[data["temp","population","numbers"][[50;;150,;;300]],AspectRatio1/GoldenRatio,ImageSizeLarge]

Out[]=

Let’s extract the population data for the US.

GeoBoundingBox

Continental US states

ADMINISTRATIVE DIVISIONS



Out[]=

{GeoPosition[{24.5455,-124.733}],GeoPosition[{49.3844,-66.9499}]}

The following equations convert the grid to latitude and longitude.

In[]:=

lon[n_]:=.25-180+(n-1)30/60lat[n_]:=-.25+90-(n-1)30/60

{n/.Solve[lon[n]-124.7329849,n,Reals],n/.Solve[lon[n]-66.9499,n,Reals],n/.Solve[lat[n]25.5455,n,Reals],n/.Solve[lat[n]49.3844,n,Reals]}

Out[]=

{{111.034},{226.6},{129.409},{81.7312}}

Now that we now the limits in the grid let’s proceed to extract the data of interest.

In[]:=

data["temp","population","assocGeoPosition"]=Association[DeleteCases[Flatten[Table[GeoPosition[{lat[i],lon[j]}]data["temp","population","numbers"]〚i,j〛,{i,82,129},{j,111,227}],1],Alternatives[Rule[_,0.],Rule[_,0]]]];

We must remove the points outside the US polygon/

In[]:=

poly=

Continental US states

ADMINISTRATIVE DIVISIONS

["Polygon"]/.GeoPositionSequence;

points=Keys@data["temp","population","assocGeoPosition"];validPoints=Fold[Or,False,#]&/@(Transpose[Table[With[{values=points/.GeoPositionSequence},RegionMember[poly[[n]],values]],{n,Length@poly}]]);Length@validPoints

Out[]=

4231

In[]:=

data["temp","population","assocGeoPosition"]=Association[Extract[Normal[data["temp","population","assocGeoPosition"]],Position[validPoints,True]]];

In[]:=

Clear[points,validPoints]

GeoGraphicsPolygon

Continental US states

ADMINISTRATIVE DIVISIONS

,PointSize[0.0001],Point[Keys@data["temp","population","assocGeoPosition"]]

Out[]=

data["temp","population","assocGeoPosition"][[;;3]]

Out[]=

GeoPosition[{48.75,-122.25}]187237.,GeoPosition[{48.75,-121.75}]5404,GeoPosition[{48.75,-121.25}]730.653

As a side note, we can observe the distribution of the population across the different cells of the grid.

Histogram[Log[2,Values[data["temp","population","assocGeoPosition"]]],PlotTheme"Detailed",FrameTicks{{Automatic,None},{{#,2^#}&/@Range[0,20,5],None}}]

Out[]=

GeoRegionValuePlotLog[2,data["temp","population","assocGeoPosition"]],ColorFunctionScalingFalse,ColorFunctionWhich#≤2,

,#≤4,

,#≤6,

,#≤12,

,#≤16,

,#≤18,

,#≤20,LightYellow,True,LightPink&

Out[]=

The mobility layers

The tessellation procedure will divide the US surface in a metapopulation structure, each with its own subpopulations. The spatio-temporal patterns of disease spreading will be linked to the mobility flows of people between the different subpopulations. These flows will be represented by two different mobility networks.

Air Traffic Network

The Bureau of Transportation Statistics keeps air traffic statistics within the US

. Table T-100 Market (all carriers) collects the number of passengers enplaned between two airport within a period. We proceeded to download the data for the year 2018.

data["sources","airTrafficData"]=Association["home"Hyperlink["https://www.bts.gov/topics/airlines-and-airports-0"],"tables"Hyperlink["https://www.transtats.bts.gov/Tables.asp?DB_ID=111&DB_Name=Air%20Carrier%20Statistics%20%28Form%2041%20Traffic%29-%20All%20Carriers&DB_Short_Name=Air%20Carriers"]]

Out[]=

home

https://www.bts.gov/topics/airlines-and-airports-0

,tables

https://www.transtats.bts.gov/Tables.asp?D

B_ID=111&DB_Name=Air%20Carrier%20Statistics%20%28Form%2041%20Traffic%29-%20All%20Carriers&DB_Short_Name=

Air%20Carriers



In[]:=

data["rawData","t100market"]=Import[FileNameJoin[{NotebookDirectory[],"datasets","airTraffic","data","351722122_T_T100D_MARKET_ALL_CARRIER.CSV"}]];

In[]:=

data["temp","airTraffic"]=Association[];

data["temp","airTraffic","header"]=MapIndexed[{First@#2,#1}&,data["rawData","t100market"][[1]]]

Out[]=

{{1,PASSENGERS},{2,ORIGIN_AIRPORT_ID},{3,ORIGIN},{4,DEST_AIRPORT_ID},{5,DEST},{6,YEAR},{7,MONTH},{8,}}

In[]:=

data["temp","airTraffic","data"]=Rest@data["rawData","t100market"];

First step is always to take a look at the data. There is always some issue that can be found. In the following steps we’ll proceed to clean up data.

Example, remove flights with no passengers.

In[]:=

data["temp","airTraffic","data"]=DeleteCases[data["temp","airTraffic","data"],{0.,__}];

Remove Flights to and from the same airport and last field.

In[]:=

data["temp","airTraffic","data"]=Most/@DeleteCases[data["temp","airTraffic","data"],{_,_,a_,_,a_,__}];

Obtain the list of airports.

data["temp","airTraffic","data"][[;;7]]

Out[]=

{{1.,10304,ANI,12719,KLG,2018,5},{1.,10304,ANI,13196,MCG,2018,5},{1.,10313,ANV,10304,ANI,2018,5},{1.,10313,ANV,12676,KGX,2018,5},{1.,10551,BET,12709,KKI,2018,5},{1.,10551,BET,14130,PKA,2018,5},{1.,10551,BET,15862,WTL,2018,5}}

Filter only to airports within the 48 contiguous states and listed in the air traffic database.

In[]:=

data["temp","airTraffic","airports"]=Association[];

Airports within the USA.

EntityPrefetch["Airport"]

Out[]=

Success

✓

Message: Prefetch successful.

Values: 1161056



In[]:=

data["temp","airTraffic","airports","usAirports"]=EntityListFilteredEntityClass"Airport",EntityFunctionc,c

country

===

United States

COUNTRY

;

Filter out the airports with the 48 lower states.

In[]:=

data["temp","airTraffic","airports","airportCoordinates"]=EntityValue[data["temp","airTraffic","airports","usAirports"],{"Entity","Coordinates"}];

In[]:=

y=Fold[Or,False,#]&/@(Transpose[Table[With[{points=data["temp","airTraffic","airports","airportCoordinates"][[All,2]]},RegionMember[poly[[n]],points]],{n,Length@poly}]]);

In[]:=

data["temp","airTraffic","airports","usAirports"]=Part[data["temp","airTraffic","airports","airportCoordinates"][[All,1]],Flatten@Position[y,True]];

In[]:=

Clear[x,y,z]

data["temp","airTraffic","airports","usAirports"][[;;3]]

Out[]=



Goldstone (Gts) Airport

Clear Creek Ranch Airport

Northern Lite Airport



In[]:=

data["temp","airTraffic","airports","assocCodes2Entities"]=Association[DeleteCases[Rule[#1,#2]&@@@EntityValue[data["temp","airTraffic","airports","usAirports"],{"IATACode","Entity"}],Rule[_Missing,_],Infinity]];

data["temp","airTraffic","airports","assocCodes2Entities"][[;;3]]

Out[]=

CUX

Cuddihy Field Airport

,CSE

Buckhorn Ranch Airport

,JCY

LBJ Ranch Airport



Filter out the airports with the 48 lower states.

In[]:=

data["temp","airTraffic","airports","airportTraffic"]=ReverseSort[Association[DeleteCases[Normal[Query[All,Total][GroupBy[data["temp","airTraffic","data"][[All,{3,1}]],FirstLast]]]/.data["temp","airTraffic","airports","assocCodes2Entities"],Rule[_String,_]]]];

In[]:=

data["temp","airTraffic","airports","dsAirports"]=Dataset[With[{traffic=AssociationThread[{"Airport","Traffic"},#]&/@({#1,#2}&@@@Normal[data["temp","airTraffic","airports","airportTraffic"]]),coords=AssociationThread[{"Airport","Position","Coordinates","IATACode"},#]&/@({#1,#2,#3,#4}&@@@EntityValue[data["temp","airTraffic","airports","usAirports"],{"Entity","Position","Coordinates","IATACode"}])},JoinAcross[traffic,coords,"Airport"]]];

data["temp","airTraffic","airports","dsAirports"]=Query[ReverseSortBy[#"Traffic"&]][data["temp","airTraffic","airports","dsAirports"]]

Out[]=

We can observe that there are several airports that service a big metropolitan area.

GeoGraphics{PointSize[0.01],Point[Keys@data["temp","population","assocGeoPosition"]],Red,PointSize[0.02],Point[Normal[data["temp","airTraffic","airports","dsAirports"][All,"Airport"]]]},GeoCenter

New York City

CITY

,GeoRangeQuantity[200,"km"]

Out[]=

We need to take into account that big urban areas are serviced by multiple airports. We’ll aggregate airports that are within a distance of 100 km of each other. The following formula allows to quickly calculate distances between locations.

geoDist=Compile[{{pos1,_Real,1},{pos2,_Real,1}},ArcCos[Sin[pos1[[1]]*Pi/180]*Sin[pos2[[1]]*Pi/180]+Cos[pos1[[1]]*Pi/180]*Cos[pos2[[1]]*Pi/180]*Cos[pos2[[2]]*Pi/180-pos1[[2]]*Pi/180]]*6371]

Out[]=

CompiledFunction

Argument count: 2

Argument types:

{{_Real,1},{_Real,1}}



Test function

geoDist[{40.6925`,-74.168667`},{40.639751`,-73.778925`}]

Out[]=

33.3914

Calculate Distance matrix

In[]:=

x=Normal[data["temp","airTraffic","airports","dsAirports"][All,"Coordinates"]];

In[]:=

data["temp","airTraffic","airports","distanceMatrix"]=Flatten[With[{x=Normal[data["temp","airTraffic","airports","dsAirports"][All,"Coordinates"]]},Table[{i,j}->geoDist[x[[i]],x[[j]]],{i,Length@x-1},{j,i+1,Length@x}]],1];

Filter airports that are more than 100km apart.

In[]:=

data["temp","airTraffic","airports","nearbyPairs"]=data["temp","airTraffic","airports","distanceMatrix"][[Flatten[Position[data["temp","airTraffic","airports","distanceMatrix"][[All,2]],a_/;a<100]]]];

Proceed to group airports

In[]:=

y=Normal[data["temp","airTraffic","airports","dsAirports"][All,"Airport"]];results={};index=Union[Flatten[data["temp","airTraffic","airports","nearbyPairs"][[All,1]]]];Monitor[While[Length@index>0,vals={z=Intersection[index,Union[Flatten[Cases[data["temp","airTraffic","airports","nearbyPairs"][[All,1]],Alternatives[{First@index,_},{_,First@index}]]]]];If[Length@z0,Null,group=y[[z]];First@groupgroup]};results=Join[results,vals];index=Complement[index,Flatten[{z,First@index}]];],First@index]results=DeleteCases[results,Null];results=Join[results,Rule[#,{#}]&/@Complement[y,Union@Flatten@Values[results]]];data["temp","airTraffic","airports","airportGroups"]=Association[results];Clear[x,y,vals,results,index]

(*GeoGraphics[{PointSize@Large,Red,Tooltip[Point[#],#]&/@Flatten@Keys[data["temp","airTraffic","airports","airportGroups"]]},GeoScaleBarPlaced[{"Metric","Imperial"},{Right,Bottom}],ImageSizeLarge]*)

In[]:=

$UnitSystem="Metric";

Map the airport hubs

GeoGraphics[{PointSize[0.001],Point[Keys@data["temp","population","assocGeoPosition"]],Red,PointSize[0.005],Tooltip[Point[#],#]&/@Keys@data["temp","airTraffic","airports","airportGroups"]},ImageSize800]

Out[]=

We will now aggregate all traffic data to the corresponding airport hubs. Let’s work with IATA codes instead of entities.

In[]:=

data["temp","airTraffic","airports","assocEntities2Codes"]=Association[Rule[#2,#1]&@@@Normal@data["temp","airTraffic","airports","assocCodes2Entities"]];

In[]:=

data["temp","airTraffic","airports","IATACodeGroupMembership"]=Association[Normal[data["temp","airTraffic","airports","airportGroups"]]/.data["temp","airTraffic","airports","assocEntities2Codes"]];

data["temp","airTraffic","airports","IATACodeGroupMembership"][[;;3]]

Out[]=

ATL{ATL,CSG,MCN,PDK,AHN,ANB,MGE,FTY,GVL,AUO},ORD{ORD,MDW,MKE,SBN,RFD,GYY,PWK,DPA,UGN,LOT,ENW,BEH,UES},LAX{LAX,SNA,BUR,ONT,LGB,SBA,RIV,CLD,HHR,VNY,VCV,SMO,OXR,EDW,TOA}

The following association points to the associated airport hub for each airport

In[]:=

data["temp","airTraffic","airports","IATACode2HubCode"]=Association@Flatten[Thread[Rule[#2,#1]]&@@@Normal[data["temp","airTraffic","airports","IATACodeGroupMembership"]]];

Replace airports for hubs in the dataset. Let’s first leave only the airports that we have identified.

In[]:=

x=DeleteCases[(data["temp","airTraffic","data"][[All,{3,5,1}]])/.data["temp","airTraffic","airports","assocCodes2Entities"],Alternatives[{_String,_,_},{_,_String,_}]]/.data["temp","airTraffic","airports","assocEntities2Codes"]/.data["temp","airTraffic","airports","IATACode2HubCode"];

In[]:=

data["temp","airTraffic","airports","IATACodeGroupIndex"]=Association[MapIndexed[Rule[#1,First@#2]&,data["temp","airTraffic","airports","IATACodeGroupMembership"]//Keys]];

data["temp","airTraffic","airports","airTrafficHubMatrix"]=SparseArray[DeleteCases[SortBy[Normal[Query[GroupBy[{#[[1]],#[[2]]}&],Total,#[[3]]&][x]]/.data["temp","airTraffic","airports","IATACodeGroupIndex"],First],Alternatives[{a_,a_}_,{_String,_}_,{_,_String}_],Infinity]]

Out[]=

SparseArray

Specified elements: 5680

Dimensions: {167,167}



data["temp","airTraffic","airports"]//Keys

Out[]=

{usAirports,airportCoordinates,assocCodes2Entities,airportTraffic,dsAirports,distanceMatrix,nearbyPairs,airportGroups,assocEntities2Codes,IATACodeGroupMembership,IATACode2HubCode,IATACodeGroupIndex,airTrafficHubMatrix,airportHubCoordinates,airportHubIndex}

In[]:=

Clear[x]

Checking that the air traffic mobility network looks correct.

ShowGraphicsFaceForm[LightGray],EdgeForm[Brown],

Continental US states