[R-sig-Geo] Units of centroid point coordinates generated byget.Pcent and calcCentroid?
Rick Reeves
reeves at nceas.ucsb.edu
Wed Aug 29 19:41:57 CEST 2007
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:
> 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
>>
>
--
Rick Reeves
Scientific Programmer / Analyst
National Center for Ecological Analysis and Synthesis
UC Santa Barbara
reeves at nceas.ucsb.edu
www.nceas.ucsb.edu
805 892 2533
-------------- next part --------------
A non-text attachment was scrubbed...
Name: reeves.vcf
Type: text/x-vcard
Size: 339 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20070829/e6851ca2/attachment.vcf>
More information about the R-sig-Geo
mailing list