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 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.

share|improve this question
1  
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

3 Answers 3

up vote 8 down vote accepted

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

Mathematica graphics

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"]

Mathematica graphics

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.

Mathematica graphics

The solution is to map Reverse on these coordinates:

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

Mathematica graphics

share|improve this answer
    
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}]]]

output

share|improve this answer
1  
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

Lets start with a map of the USA:

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

enter image description here

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]]]]

enter image description here

share|improve this answer
    
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

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.