[R-sig-Geo] Units of centroid point coordinates generated byget.Pcent and calcCentroid?

[Roger Bivand]

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

-- 
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: Roger.Bivand at nhh.no