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

[Michael Fuller]

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