6
$\begingroup$

I'm importing the NGA ASAM piracy database into Mathematica for analysis.

asamjson=Import["http://msi.nga.mil/MSI_JWS/ASAM_JSON/getJSON?typename=\
DateRange_AllRefNumbers&fromDate=19900101&toDate=20140801", 
 "JSON"];
Needs["GeneralUtilities`"]
asamdataset = Dataset[ToAssociations@asamjson];
piracyLocations = 
 Map[GeoPosition[ ToExpression@{#["lat"], #["lng"] }]   &, 
  asamdataset];

I can import the data fine, and get lat/longs (for some reason they're included twice in the ASAM data). I can then draw them with GeoListPlot and get map of all the piracy events.

GeoListPlot[Take[piracyLocations, 300]]

However, when I then go to make a heat map histogram things fall apart. What displayed nicely now shows up with a wildly out of skew map. It looks like I need to choose some different form of projection, but which?

dat = Map[ ToExpression@{#["lat"], #["lng"]}&,Take[asamdataset, 4000]];

histo = SmoothDensityHistogram[dat, .3, 
 ColorFunction -> "SunsetColors", PlotPoints -> 150, Mesh -> 20, 
 MeshStyle -> Opacity[.1], PlotRange -> All];

Show[histo, 
 Graphics[{FaceForm[], EdgeForm[Directive[White, Opacity[.7]]], 
 CountryData["World", "Polygon"]}]]

How do I project my heatmap of piracy events correctly onto the world?

EDIT: My first problem is that Geoposition and related functions take Lat/Long pairs. Which is Y/X, whereas other Mathematica functions take X/Y pairs.

$\endgroup$
1
  • 2
    $\begingroup$ The reason for not being able to plot more than 4000 points is that there are two null s in the piracyLocations. In[1079]:= Position[Normal@piracyLocations,null] Out[1079]= {{4335,1,2},{4391,1,1}} These cause the plot to freeze. $\endgroup$
    – user22228
    Nov 18, 2014 at 18:48

3 Answers 3

6
$\begingroup$

You can also exploit the GeoGraphics in version 10. Using some of your variables (noting there seems to be an issue with plotting more than 4000 points or generating density histogrm with more points...that I do not fully understand):

Getting a feel:

GeoGraphics[{Red, PointSize[0.005], 
  Point[(Normal@piracyLocations)[[1 ;; 4000]]]}]

enter image description here

Now generating heat map:

sdh = SmoothDensityHistogram[
   Reverse[First[#]] & /@ (Normal@piracyLocations)[[1 ;; 4000]], 
   ColorFunction -> "SunsetColors", 
   PlotRange -> {{-180, 180}, {-90, 90}}];
GeoGraphics[{Opacity[0.6], sdh[[1]]}, GeoBackground -> "ReliefMap"]

enter image description here

$\endgroup$
1
  • $\begingroup$ Why do does the histogram smooth so much on the X axis? Even if I use a PDF smoother it seems wildly inaccurate. I'm sure it isn't, but something seems off there. $\endgroup$
    – user4860
    Aug 4, 2014 at 0:44
8
$\begingroup$

You have mixed up latitudes and longitudes. So first of all change ToExpression@{#["lat"], #["lng"] } into ToExpression@{#["lng"], #["lat"] }. The second problem you have is that when you combine graphics with Show all the common options are taken to be that of the first graphics which in your case is the histogram, but what you really want is the aspect ratio, size etc. of the map. Just to demonstrate I did a quick fix:

Show[Graphics[{FaceForm[], EdgeForm[Directive[White, Opacity[.7]]], 
   CountryData["World", "Polygon"]}], histo, 
 Graphics[{FaceForm[], EdgeForm[Directive[White, Opacity[.7]]], 
   CountryData["World", "Polygon"]}]]

Heat map

$\endgroup$
0
4
$\begingroup$

Off Topic, but to show another approach on showing the concentration of events using Histogram3D, which can convey the concentration of events a wee bit better.

asamjson = 
  Import["http://msi.nga.mil/MSI_JWS/ASAM_JSON/getJSON?typename=\
DateRange_AllRefNumbers&fromDate=19900101&toDate=20140801", "JSON"];
Needs["GeneralUtilities`"]
asamdataset = Dataset[ToAssociations@asamjson];
coord10 = 
  Normal@Map[ToExpression@{#["lat"], #["lng"]} &, asamdataset];
coord10[[All, 2]] = -coord10[[All, 2]];
geoRange = GeoRange -> {{-45, 68}, {-180, 180}};
background = 
  Graphics3D[{Texture[GeoGraphics[geoRange]], EdgeForm[], 
    Polygon[{{-45., 180., 0.}, {-45., -180., 0.}, {68., -180., 
       0.}, {68., 180., 0.}}, 
     VertexTextureCoordinates -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}}]}, 
   Lighting -> "Neutral"];
Show[Histogram3D[coord10, {1}, BoxRatios -> {68 + 45, 360, 90}, 
  Axes -> False, ChartStyle -> Red, Boxed -> False, 
  ViewPoint -> {-1, 0, 3}, ImageSize -> 800], background]

Mathematica graphics

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.