[R-sig-Geo] Problem with aggregating UScensus2000 FL data

[Roger Bivand]

On Mon, 1 Aug 2011, Roger Bivand wrote:

> On Mon, 1 Aug 2011, Ista Zahn wrote:
>
>> Hi Roger,
>> 
>> On Mon, Aug 1, 2011 at 2:38 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>>> On Mon, 1 Aug 2011, Ista Zahn wrote:
>>> 
>>>> Oops, I got the examples mixed up. The Alabama example works, the
>>>> Florida example does not. I also should have given the version info,
>>>> as this may well be version specific.
>>> 
>>> You also need to provide the link to the data file, and the code prior to
>>> what you show. The function county() is not defined anywhere, consequently
>>> the attribute fields are nothing like tract boundaries such as those on:
>> 
>> Sorry I wasn't clear. The data and the county function are in the
>> UScensus2000 packages. Let me try again. A complete example leading to
>> the error is
>> 
>> library(UScensus2000)
>
> Perhaps contact the package maintainer? It is very possible that the package 
> ships with unclean topologies, isn't it? I'll try to look later.

library(UScensus2000)
library(maptools)
data(countyfips)
tmp <- county(fips = gsub("^..", "", countyfips$fips[countyfips$acronym == 
"fl"]), state = "fl", level = "tract")
dat <- as.data.frame(tmp)
tmp.county <- unionSpatialPolygons(tmp, dat$county)
# fails on (-82.6787 28.4335) == (-82.6787 28.4335)
getScale()
setScale(1e+09)
tmp.county <- unionSpatialPolygons(tmp, dat$county)
# works

The issue is that GEOS assumes that it is working on planar coordinates, 
and allows plenty of trailing significant digits before truncating. With 
geographical coordinates, the nth trailing digit may be the one that 
disambiguates two points. Increasing the scale by one order of magnitude 
is sufficient in this case. I do not see the coordinates you stated 
initially.

Roger

>
> Roger
>
>> library(maptools)
>> data(countyfips)
>> 
>> tmp <- county(fips = gsub("^..", "",
>> countyfips$fips[countyfips$acronym == "fl"]), state = "fl", level =
>> "tract")
>> dat <- as.data.frame(tmp)
>> tmp.county <- unionSpatialPolygons(tmp, dat$county)
>> 
>>> 
>>> http://www.census.gov/geo/www/cob/tr2000.html
>>> 
>>> The lines do intersect (touch), and should not in a topologically correct
>>> setting; the coordinates seem odd - what is the projection? Did you run
>>> checkPolygonsHoles() first?
>> 
>> I would if I could figure out how to convert a SpatialPolygons object
>> to a Polygons object... I'm just getting started with sp.
>> 
>> Thanks,
>> Ista
>> 
>>> 
>>> Roger
>>> 
>>>> 
>>>> R version 2.13.1 (2011-07-08)
>>>> Platform: i386-pc-mingw32/i386 (32-bit)
>>>> 
>>>> locale:
>>>> [1] LC_COLLATE=English_United States.1252
>>>> [2] LC_CTYPE=English_United States.1252
>>>> [3] LC_MONETARY=English_United States.1252
>>>> [4] LC_NUMERIC=C
>>>> [5] LC_TIME=English_United States.1252
>>>> 
>>>> attached base packages:
>>>> [1] grid      stats     graphics  grDevices utils     datasets  methods
>>>> [8] base
>>>> 
>>>> other attached packages:
>>>> [1] rgeos_0.1-8             stringr_0.5             quadprog_1.5-4
>>>> [4] directlabels_1.3        UScensus2000_1.00       gpclib_1.5-1
>>>> [7] UScensus2000cdp_0.03    UScensus2000blkgrp_0.03 
>>>> UScensus2000tract_0.03
>>>> [10] rgdal_0.7-1             maptools_0.8-9          lattice_0.19-30
>>>> [13] sp_0.9-84               foreign_0.8-45          ggplot2_0.8.9
>>>> [16] proto_0.3-9.2           reshape_0.8.4           plyr_1.5.2
>>>> 
>>>> loaded via a namespace (and not attached):
>>>> [1] digest_0.5.0 tools_2.13.1
>>>>> 
>>>> 
>>>> On Mon, Aug 1, 2011 at 11:09 AM, Ista Zahn <izahn at psych.rochester.edu>
>>>> wrote:
>>>>> 
>>>>> Hi all,
>>>>> I am trying to aggregate the UScensus data for Florida up to the zip
>>>>> code level using unionSpatialPolygons from the maptools package. This
>>>>> works for other states, but not for Florida:
>>>>> 
>>>>> ## Works
>>>>> tmp <- county(fips = gsub("^..", "",
>>>>> countyfips$fips[countyfips$acronym == "fl"]), state = "fl", level =
>>>>> "tract")
>>>>> dat <- as.data.frame(tmp)
>>>>> tmp.county <- unionSpatialPolygons(tmp, dat$county)
>>>>> 
>>>>> ## Returns error
>>>>> tmp <- county(fips = gsub("^..", "",
>>>>> countyfips$fips[countyfips$acronym == "al"]), state = "al", level =
>>>>> "tract")
>>>>> dat <- as.data.frame(tmp)
>>>>> tmp.county <- unionSpatialPolygons(tmp, dat$county)
>>>>> Error in TopologyFunc(groupID(spgeom[ids[[i]]], id[ids[[i]]]),
>>>>> names(ids)[i],  :
>>>>>  TopologyException: found non-noded intersection between LINESTRING
>>>>> (-0.640183 0.433522, -0.678743 0.433456) and LINESTRING (-0.678743
>>>>> 0.433456, -0.677839 0.434367) at -0.678743 0.433456
>>>>> 
>>>>> Does anyone have a suggestion about how to do this?
>>>>> 
>>>>> Thanks!
>>>>> 
>>>>> Ista
>>>>> --
>>>>> Ista Zahn
>>>>> Graduate student
>>>>> University of Rochester
>>>>> Department of Clinical and Social Psychology
>>>>> http://yourpsyche.org
>>>>> 
>>>> 
>>>> 
>>>> 
>>>> 
>>> 
>>> --
>>> Roger Bivand
>>> Department of Economics, NHH Norwegian School of Economics,
>>> 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
>>> 
>> 
>> 
>> 
>> 
>
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
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