[R-sig-Geo] Units of centroid point coordinates generated byget.Pcent and calcCentroid?
Roger Bivand
Roger.Bivand at nhh.no
Wed Aug 29 20:20:16 CEST 2007
On Wed, 29 Aug 2007, Rick Reeves wrote:
> 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.
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
>
> Regards,
> RR
>
>
>
>
>
>
>
>
>
> Roger Bivand wrote:
>> 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
More information about the R-sig-Geo
mailing list