Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question

3 Answers 3

up vote 3 down vote accepted

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

share|improve this answer
    
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. –  Steven Aug 4 at 0:44

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

share|improve this answer

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.