Skip to content

Units of centroid point coordinates generated by get.Pcent and calcCentroid?

5 messages · Rick Reeves, Michael Sumner, Roger Bivand

Original
Rick Reeves
# 
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
Michael Sumner
# 
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.

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:
Roger Bivand
#
On Wed, 29 Aug 2007, Michael Sumner wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Rick Reeves
    

    

      
      #
    

  


  

      
Thank you both for your answers:

Here is what I think that you are trying to convey:

1) The centroid points do not match the map projection because the 
projection attributes for the shapefile are poorly defined.
2) If I correctly define the shape file's projection attributes, then 
both get.Pcent() and calcCentroid() generate centroids
    in units that match the projection of the shape file.
 
One additional fact: this shape file does have a .prj file:

PROJCS["North_America_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",20.0],
PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",40.0],UNIT["Meter",1.0]]

I made the assumption that:
MapPolysMapClass = read.shape("FiveNWStateCounties")

...would read the .prj file. I guess not, from reading your message.

Do the projection settings specified by spTransform() override the 
contents of the .prj file (when it exists)? I think the answer is yes.

I will try 1) using the OGR methods, 2) Using the spTransform methods to 
explicitly define the projection parameters.

Regards,
RR

        
Roger Bivand wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 29 Aug 2007, Rick Reeves wrote:

        

            

      
      
      
    

        
You are right that none of the functions in maptools read the *.prj file. 
The readOGR() function in rgdal does read it, so if that suits you, no 
transformation is required. However, if you want to check the centroids in 
the file, use spTransform() to go from the projected CRS read by readOGR() 
to geographical coordinates. If the fine details matter, you'd need to 
know which ellipsoid and datum were used in the polygons used to calculate 
the centroid values you have.

Roger