# plotting a {{Latitude,Longitude}} list on top of a GIS shapefile

I have a dataset composed of points specified in {Latitude,Longitude} format. Call it latlong for the purposes of this question. It displays properly with

Graphics[{Point /@
Map[GeoGridPosition[GeoPosition[#],
"Mercator"][[1]] &, {latlong}, {2}]}]


I also have a shape file that contains the boundaries of the counties in which the latlong points reside. If I query Import[demo.shp", "CoordinateSystemInformation"], I get

"GEOGCS" -> {"GCS_North_American_1983",
"DATUM" -> {"North_American_Datum_1983",
"SPHEROID" -> {"GRS_1980", 6.37814*10^6, 298.257}},
"PRIMEM" -> {"Greenwich", 0.}, "UNIT" -> {"Degree", 0.0174533}}


and if I just Import the shape file, it displays correctly on screen.

I would like to display the points from latlong within the county borders from demo.shp. I have not been able to figure out how.

I appreciate that this is as much a GIS question as it is a Mathematica question, but I hope somebody here can help.

-
Can you upload demo.shp somewhere? Or just include a small sample of it –  ssch Jan 23 '13 at 19:22
Assuming your latlong dataset and your SHP use the same datum I don't see a problem in combining the graphics resulting from both. What precisely is your problem? –  Sjoerd C. de Vries Jan 23 '13 at 19:50
The map invariably appears as a dot in one corner of the graph, and the points as a dot in the other corner. I suspect I'm messing up the projections. The data can be found at box.com/s/qw69but73asus7mig3j3 Let me know if you have any trouble accessing that folder, and thanks in advance. –  Michael Stern Jan 23 '13 at 20:12
Mathematica complains that it isn't an SHP file. (BTW don't forget to use the @name convention if you want to ping someone). –  Sjoerd C. de Vries Jan 23 '13 at 20:21
for a valid shapefile:.shp, .shx, .dbf are three Mandatory files –  Tuku Jan 23 '13 at 20:23

The most likely cause of your problem is that you are swapping lat and long coordinates.

Your latlong definition in the .NB file provided goes like this:

latlong = {{40.660323, -73.997952}, {40.660489, -73.99822}, {40.654365, -74.004113},...


That's New York. Simply plotting without projection gives:

m1=Graphics@Point@latlong


I seem to see Long Island and some other familiar features. Looks like a subway map.

Here is a publicly available SHP map of the subway lines:

Import["https://wfs.gc.cuny.edu/SRomalewski/MTA_GISdata/June2010_update/nyctsubwayroutes_100627.zip", "SHP"]


Extracting some coordinates from that:

Cases[map, {_Real, _}, Infinity]


{{-74.015047, 40.703577}, {-74.015028, 40.703214}, {-74.014889, 40.702506}, ...

You can see that this has the latlongs in the reverse order. Plotting these points with your other latlong gives precisely the result you describe above.

The solution is to map Reverse on these coordinates:

m2 = Graphics[{Red, Point@(Reverse /@ Cases[map, {_Real, _}, Infinity])}];
Show[m1,m2]


-
The fact that the answer was so straight-forward doesn't make me any less impressed with how well you figured it out. I should have noticed that my dataset was placing Far Rockaway in the upper left, but somehow I missed it. Thanks! –  Michael Stern Jan 23 '13 at 21:02
nice detective work +1 –  Mike Honeychurch Jan 23 '13 at 21:07
bg = Import["sample.shp"]
Show[bg, Graphics[Point@Reverse[latlong, {2}]]]


-
Totally correct, and thanks. Sjoerd C. de Vries beat you by a few minutes so I gave him the points. –  Michael Stern Jan 23 '13 at 21:03

 plot1 = Import["http://exampledata.wolfram.com/usamap.zip","Graphics"]


Now we can add points on top of the map easily:

 points = Transpose@{RandomReal[{-100, -90}, 50], RandomReal[{35, 47}, 50]};

Show[plot1,  ListPlot[points, PlotStyle -> Directive[Red, PointSize[Large]]]]


-
I might be wrong, but I assume the OP's problem is not in combining the graphics but it in taking care they use the same coordinate system. I don't believe you address that problem here. –  Sjoerd C. de Vries Jan 23 '13 at 20:00
correct, each renders correctly independently but placement and scale are very wrong when I combine them. I presume I'm messing up the coordinate system. –  Michael Stern Jan 23 '13 at 20:13