Exploring the globe
** Please note** The following example is for illustration only. For the most recent examples and documentation, please visit our GitHub page.
Suppose you want to explore world geography vis-a-vis networks of neighboring cities. A wealth of data is available online, and for this exercise we will use a catalog of 7322 cities from simplemaps.com. If you are unable to download the .csv file from the hyperlink, you can copy and paste to a text editor and save with extension .csv (alternatively, a copy can be found in the Eirene.jl folder). Identify the <filepath> to the download, and load to a Julia REPL using the Eirene ezread wrapper:
julia> include("<filepathtoEirene>")
julia> a = ezread("<filepathtocsvfile>")
7323x9 Array{Any,2}:
"city" "city_ascii" "lat" "lng" … "country" "iso2" "iso3" "province"
"Qal eh-ye Now" "Qal eh-ye" 34.983 63.1333 "Afghanistan" "AF" "AFG" "Badghis"
"Chaghcharan" "Chaghcharan" 34.5167 65.25 "Afghanistan" "AF" "AFG" "Ghor"
"Lashkar Gah" "Lashkar Gah" 31.583 64.36 "Afghanistan" "AF" "AFG" "Hilmand"
"Zaranj" "Zaranj" 31.112 61.887 "Afghanistan" "AF" "AFG" "Nimroz"
"Tarin Kowt" "Tarin Kowt" 32.6333 65.8667 … "Afghanistan" "AF" "AFG" "Uruzgan"
"Zareh Sharan" "Zareh Sharan" 32.85 68.4167 "Afghanistan" "AF" "AFG" "Paktika"
⋮ ⋱ ⋮
"Gweru" "Gweru" -19.45 29.82 "Zimbabwe" "ZW" "ZWE" "Midlands"
"Mutare" "Mutare" -18.97 32.65 "Zimbabwe" "ZW" "ZWE" "Manicaland"
"Kadoma" "Kadoma" -18.33 29.9099 "Zimbabwe" "ZW" "ZWE" "Mashonaland West"
"Chitungwiza" "Chitungwiza" -18.0 31.1 … "Zimbabwe" "ZW" "ZWE" "Harare"
"Harare" "Harare" -17.8178 31.0447 "Zimbabwe" "ZW" "ZWE" "Harare"
"Bulawayo" "Bulawayo" -20.17 28.58 "Zimbabwe" "ZW" "ZWE" "Bulawayo"
A comment on the meaning of Array{Any,2} in the third line: Julia differentiates arrays based on the type (or class) of elements they can contain. An n-dimensional array that is formally declared to have elements of type T is called an array of type Array{T,n}. A 5x5 matrix with 64-bit floating point entries can be realized either as an array of type Array{Float64,2} or, as the name suggests, as an array of type Array{Any,2}. Many functions in Julia are specialized to work with numerical arrays, so it's in our interest to extract the numeric part we are interested in, namely columns 3 and 4, sans headers, and convert to an array of type Array{Float64,2} before proceeding.
julia> b = a[2:end,3:4]
7322x2 Array{Any,2}:
34.983 63.1333
34.5167 65.25
31.583 64.36
⋮
-17.8178 31.0447
-20.17 28.58
julia> b = convert(Array{Float64,2},b)
7322x2 Array{Float64,2}:
34.983 63.1333
34.5167 65.25
31.583 64.36
⋮
-17.8178 31.0447
-20.17 28.58
Eirene has a built-in function to convert spherical coordinates (in degrees, with fixed radius 1) to 3D Euclidean coordinates. The rowsare keyword argument determines wether rows are treated as points or dimensions.
julia> c = latlon2euc(b,rowsare="points")
3x7322 Array{Float64,2}:
0.370265 0.344959 0.368622 0.403432 0.344318 0.309031 … 0.796253 0.822829 0.814358 0.81567 0.824296
0.730885 0.748275 0.767998 0.755149 0.768533 0.781189 0.510205 0.473338 0.491252 0.490971 0.449048
0.573333 0.566646 0.523733 0.516713 0.53926 0.542442 -0.325073 -0.31449 -0.309017 -0.305991 -0.344807
Now that we have a point cloud we can take a preliminary look at its shape and scale with the ezplot_pjs wrapper. Note you can always select "deny" when PlotlyJS asks for web access, without penalty.
julia> ezplot_pjs(c)
A cursory inspection shows that a number of interesting features appear at or below epsilon = 0.15, so we'll use that as our initial cutoff. To pass city names to the eirene function, extract column 2 of array a, sans header, and store in a 1-dimensional array of type Array{Any}. Due to potential errors resulting from nonstandard character strings, it's good practice to use the built-in label sanitzer ezlabel to clean this column before assigning to d. The wrapper will replace any element of the column that cannot be expressed as an ACIIString with the number corresponding to its row. Aside: n can be omitted from the expression Array{T,n} when n=1.
julia> d = ezlabel(a[2:end,2])
7322-element Array{Any,1}:
"Qal eh-ye"
"Chaghcharan"
"Lashkar Gah"
⋮
"Harare"
"Bulawayo"
With our inputs in order, it's time to call eirene. The calculation should take around 2.5GB of memory and 2 min to complete.
julia> C = eirene(c,rowsare = "dimensions",upperlim = 0.15,pointlabels=d)
elapsed time: 118.869523387 seconds
Dict{ASCIIString,Any} with 14 entries:
"symmat" => 7322x7322 Array{Int64,2}:…
"filtration" => Any[[515055,515055,515055,515055,515055,515055,515055,515055,515055,515055 … 515055,51505…
"lowlab0" => Any[Int64[],[1,2,3,4,5,6,7,8,9,10 … 7313,7314,7315,7316,7317,7318,7319,7320,7321,7322],[3…
"firstv" => Any[[1,2,3,4,5,6,7,8,9,10 … 7314,7315,7316,7317,7318,7319,7320,7321,7322,7323],[1,346,688…
"filtrationtranslator" => [0.149999,0.149999,0.149999,0.149999,0.149999,0.149999,0.149999,0.149998,0.149998,0.149998 …
⋮ => ⋮
The packge JLD can be used to save this output. Run Pkg.add("JLD") if you have not done so already, and enter
julia> JLD.save("demography.jld","C",C)
The file will be saved to your current working directory. To recover what you have saved, run
julia> X = JLD.load("demography.jld")
This returns a dictionary object X with one key, "C". Running
julia> C = X["C"]
will leave C unchanged.
We did not specify bettimax, so only the persistence modules in dimensions 0 and 1 are computed. To view the persistence diagram in dimension 1, enter
julia> plotpersistencediagram_pjs(C)
Hovering over a point in the diagram will display the identification number of the corresponding persistent homology class, together with the size the cycle representative Eirene computed for it. Fix a feature of interest - say number 1757 - and plot the corresponding cycle, noting that only vertices will appear in the figure (cells of dimension 1 and higher are generally too numerous to plot efficiently).
Note: The cost of displaying text labels for individual points in the Plotly platform is higher than that of displaying points alone. Consider using showlabels=false, if your graphics card has any difficulties.
julia> plotclassrep_pjs(C,class = 1757)
Class 1757, shown above, is among my favorites: the Himalayan branch of the silk road is clearly visible on its South West arc; to the North it follows the Trans-Siberian railroad from Moscow in the west to Vladivostok on the Sea of Japan, by way of Omsk, Irktsk, and Chita; from there, it follows the connecting route from Beijing to Hong Kong, passing Nanning, Hanoi, and close to Ho Chi Min on its way to Bangkok (not all of these appear in the cycle itself, but they are easily spotted nearby). With a little exploration you can find large features formed by the Sahara Desert, Tapajos River, Falkland Islands, South China Sea, Guam, and Hudson Bay, and smaller ones shaped by the Andes mountains and Gulf of Mexico - all through the proxy of urban development.
To plot the same class-representative without text labels, use
julia> plotclassrep_pjs(C,class = 1757,showlabels=false)
To plot only the cities incident to the cycle representative, use
julia> plotclassrep_pjs(C,class = 1757,showcloud=false)
To plot these using MDS, use
julia> plotclassrep_pjs(C,class = 1757,showcloud=false,coords="mds",embeddingobj="hop")
To make any further changes to the appearance of your plot, share it with others, or export to still-frames, you can upload to the interactive, open-source Plotly web API. This requires the setup of a personal account beforehand (see the github documentation for Plotly.jl for instructions), but once this has been done and the connection between your Julia REPL and Plotly account has been established, uploading is quite simple: create an object p for the PlotlyJS figure, and call Plotly.post
julia> p = plotclassrep_pjs(C,class = 1757)
julia> Plotly.post(p)
Plotly.RemotePlot(URI(https://plot.ly/~henselmonster/78))
The figure will appear under the My Charts tab of your plotly library.
We hope to add more examples like this one in the future, but there is much left to learn from the Simplemaps data set. It is instructive, for instance, to observe which features change, and how, under subsampling. If you have a fun, interesting, or educational example you'd like to share with others, please let us know!
