[R-sig-Geo] geo-newbie Q: overlay lat-long data onto projected map?

Tue Jul 18 17:22:17 CEST 2006

Michael, on http://r-spatial.sourceforge.net you'll find a package (and 
installation instructions) called spmaps, which converts an object from 
the maps package (and data base) to an sp object. If needed at all...
--
Edzer

Michael Fuller wrote:

>This is awesome help! Thank you very much Tom!
>I will give it a go.
>Mike
>
>On Jul 17, 2006, at 10:04 PM, Mulholland, Tom wrote:
>
>  
>
>>Yes you are correct in your assumption that if you are not using  
>>the "maps" package you need to have your own boundary files. There  
>>is the "National Atlas of the United States" produced by the USGS  
>>which I recall has downloads of a variety of shapefiles. http://www- 
>>atlas.usgs.gov/ is the mainpage but administrative boundaries are  
>>on http://www-atlas.usgs.gov/atlasftp.html? 
>>openChapters=chpbound#chpbound in compressed files. I'm sure I have  
>>previously found them uncompressed, but it's not too hard to  
>>decompress tham as you only have to do it once.
>>
>>
>>So once you have the shapefiles you can read them using the  
>>maptools package with the readShapePoly set of commands.
>>
>>such as
>>
>>US <- readShapePoly("name of shape file")
>>#make sure you have the correct transformation
>># set the CRS for the shapefile you have downloaded.
>># You need to check the metadata associated with your shapefile
>>
>># I typically use sommething like
>>#    proj4string(x) <- CRS("+proj=longlat +datum=NAD83")
>># Then reproject it into the format you wnat to use
>>#    x <- transform(x,CRS("+proj=utm +zone=50"))
>>#    }
>>
>># Having done that
>>
>>plot(US,...)
>>
>>#Then you add your data to the map
>>plot(S, pch = 20, add = TRUE)
>>
>>Yes you can extract the coordinates from maps
>>
>>    
>>
>>>test <- map("world","UK",xlim = c(-8.08932,3.4), ylim = c 
>>>(49.11741,61.54963))
>>>str(test)
>>>      
>>>
>>List of 4
>> $ x    : num [1:296] -6.24 -6.30 -6.23 -6.38 -6.46 ...
>> $ y    : num [1:296] 58.5 58.3 58.2 58.2 58.1 ...
>> $ range: num [1:4] -8.09  1.73 49.25 60.84
>> $ names: chr [1:25] "UK:Scotland:Isle of Lewis" "UK:Guernsey"  
>>"UK:Great Britain" "UK:Scotland:Shetland Islands:Unst" ...
>> - attr(*, "class")= chr "map"
>>
>>but it looks like you would them have to process them into an sp  
>>format. However the world map is not particularly detailed. It's  
>>possible the US county data is of more precision.
>>
>>If you just use maps then you should be able to use "points" to  
>>plot your data
>>
>>map('usa')
>>points(S$x, S$y, pch = 20)
>>
>>Not having delved into the map package I don't know exactly how it  
>>handles projections but it looks like
>>
>>map('usa', projection = "yourCRS")
>>
>>would be sufficient.
>>
>>I hope this helps.
>>
>>Tom Mulholland
>>Senior Demographer
>>Applied Research and Modelling
>>State and Regional Policy
>>Department for Planning and Infrastructure
>>Perth, Western Australia
>>+61 (08) 9264 7936
>>
>>
>>    
>>
>>>-----Original Message-----
>>>From: r-sig-geo-bounces at stat.math.ethz.ch
>>>[mailto:r-sig-geo-bounces at stat.math.ethz.ch]On Behalf Of
>>>Michael Fuller
>>>Sent: Tuesday, 18 July 2006 5:56 AM
>>>To: r-sig-geo at stat.math.ethz.ch
>>>Subject: [R-sig-Geo] geo-newbie Q: overlay lat-long data onto
>>>projected
>>>map?
>>>
>>>
>>>I've spent the last 5 days trying to figure out how to overlay the
>>>positions of geographic data on top of a map. I've scoured the man
>>>pages for the maps, mapproj, maptools, mapdata, and the huge sp
>>>package. I've experimented with countless permutations of commands.
>>>So far I've learned a lot about map projections, but I still have no
>>>ideas how to display my data the way I want! (sigh)
>>>
>>>The data are the latitude-longitude coordinates of 1,599
>>>survey sites
>>>located across the US. I want to plot these sites on a map of
>>>the US.
>>>Seems like a simple task. I've tried converting them to a
>>>SpatialPoints object and worked out how to project the sites using a
>>>given map projection (i.e. Bonne). But how do I correctly overlay
>>>(i.e. display) the data points onto a similarly projected map of the
>>>US? It might be possible to do it using overlay {sp} but for that it
>>>seems I would need the coordinates of the geopolitical boundary for
>>>the map of the US. Can these values be extracted somehow from the US
>>>map?
>>>
>>>Here's how I generated the SpatialPoints object:
>>>
>>>#x = longitude values in decimal degrees
>>>#y = latitude values in decimal degrees
>>>#create SpatialPoints object using CRS-class Bonne projection
>>>centered on longitude 100
>>>S <- SpatialPoints(cbind(x,y),CRS("+proj=bonne +lon_0=100"))
>>>
>>>#this next part plots just the data points (cloud of points)
>>>plot(S,pch=20)
>>>
>>>How do I overlay the same points on top of a map of the US such that
>>>they are correctly positioned wrt coordinates?
>>>
>>>Also, the R-sig-geo info page does show an URL for a searchable
>>>archive of R-sig-geo. Does one exist?
>>>
>>>Many thanks for your time.
>>>Mike
>>>______________________________
>>>Michael M Fuller, Ph.D.
>>>The Institute for Environmental Modeling
>>>University of Tennessee
>>>1416 Circle Drive
>>>Knoxville, TN 37996-1610
>>>
>>>
>>>
>>>
>>>	[[alternative HTML version deleted]]
>>>
>>>_______________________________________________
>>>R-sig-Geo mailing list
>>>R-sig-Geo at stat.math.ethz.ch
>>>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>      
>>>
>
>_______________________________________________
>R-sig-Geo mailing list
>R-sig-Geo at stat.math.ethz.ch
>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>  
>