On Wed, 29 Aug 2007, Michael Sumner wrote:
> Hello,
>
> In short: they are probably in UTM (metres) for whatever UTM zone your
> counties are closest to.
>
> Shape files sometimes have this information in a .prj file.
Yes, the input shapefile is not in geographical coordinates. But I don't
think it is UTM either, unless it is about 1000km north of the Equator,
which does not agree with the externally computed centroids (look like
Seattle/Vancouver?). So you'd need to establish the coordinate reference
system of the input data.
You could use the coordinates() methods for a SpatialPolygonsDataFrame
method, as read by readOGR() in rgdal, or readShapePoly() in maptools -
this returns the label point for the largest member polygon by naive area
(area in the input coordinates regardless of metrics). But to get
centroids in geographical coordinates, you can only do:
spTransform(SpatialPoints(coordinates(SPDF), proj4string=CRS(???)),
CRS("+proj=longlat +datum=WGS84"))
if you know what ??? is. You would autodetect a *.prj file by using
readOGR(), then ??? would be proj4string(SPDF), unless
is.na(proj4string(SPDF)).
Depending on the state and the date of the map, look in
EPSG <- make_EPSG() # rgdal
EPSG[grep("Washington", EPSG$note),]
or your state of choice to find a candidate. Once you have (say) added
county names as a SpatialPointsDataFrame, try to output:
writeOGR(my_points_ll, "centroids.kml", "centroids", driver="KML")
and open in Google Earth for a sanity check.
Hope this helps,
Roger
>
> If you use rgdal,
>
> library(sp)
> library(rgdal)
> d <- readOGR("FiveNWStateCounties.shp", "FiveNWStateCounties")
>
> proj4string(d)
>
> If proj4string is not NA, then you can use spTransform:
>
> d.ll <- spTransform(d, CRS("+proj=longlat"))
>
> Otherwise, you can project that matrix of coordinates with (if "m" is a
> matrix of your polygon centroids)
>
> x.ll <- project(m, "insert non-NA p4string here", inv = T)
>
> Do you have a .prj file, or other information that refers to UTM, or
> State Plane,with parameters?
>
> HTH, Mike.
>
>
> Rick Reeves wrote:
>> Hello List:
>>
>> I'm getting to know the polygon centroid - generating functions
>> get.Pcent() (maptools) and calcCentroid() (PBSmapping),
>> and would like to know the units of the centroid point coordinates.
>>
>> Working with maptools version:
>>
>> library(maptools)
>> MapPolysMapClass = read.shape("FiveNWStateCounties")
>> MapPolyCents = get.Pcent(MapPolysMapClass)
>>
>> Here is a portion of my input polygon file (CNTRDLONG/LAT computed
>> elsewhere),
>>
>> Browse[1]> MapPolysMapClass$att.data
>> STATEFIPS COUNTYFIPS COUNTYNAME LSAD LSAD_TRANS CNTRDLONG
>> CNTRDLAT F_AREA
>> 1 53 073 Whatcom 06 County
>> -121.7137 48.82591 5584761843
>> 2 53 073 Whatcom 06 County
>> -123.0569 48.98870 12378754
>> 3 30 029 Flathead 06 County
>> -114.0498 48.29519 13613065142
>> 4 30 053 Lincoln 06 County
>> -115.4052 48.54250 9518449902
>> 5 16 021 Boundary 06 County
>> -116.4629 48.76703 3309965666
>>
>> .and here are the corresponding polygon centroids
>>
>> [1,] -1762352.7 1282662.2
>> [2,] -1846616.8 1326565.8
>> [3,] -1256556.1 1099405.3
>> [4,] -1343898.6 1146606.3
>> [5,] -1410325.9 1187932.0
>>
>> The PBSmapping calcCentroid function (using PolySet inputs) generates
>> the same numbers.
>>
>> So, what units are these coordinates? I havent found the answer in the
>> maptools or PBSmapping package manuals.
>> Side question: (How) can I translate them to Lat/Long coordinates?
>>
>>
>> Thanks,
>> Rick Reeves
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo@stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: [EMAIL PROTECTED]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
