[R-sig-Geo] geo-newbie Q: overlay lat-long data onto projected map?
Michael Fuller
mmfuller at tiem.utk.edu
Tue Jul 18 17:55:31 CEST 2006
Great! Thanks Edzer
Mike
On Jul 18, 2006, at 11:22 AM, Edzer J. Pebesma wrote:
> 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
>>
>
More information about the R-sig-Geo
mailing list