[R-sig-Geo] Problem reading in polygon shapefile with readShapePoly
Steven J. Pierce
pierces1 at msu.edu
Mon Jun 22 16:15:20 CEST 2009
Roger,
I'll send you a link off-list to where you can download a zip file with
those files. In the meantime, here's a bit of context about the source of
the data: The parcel data originates from the office of a small city's tax
assessor. The PIN column in the .DBF file contains what is supposed to be a
unique ID for each parcel of land. The rest of the data largely contain the
info about the address of the property represented by the parcel, the zoning
codes and property class for the parcel, and info about the taxpayer
responsible for the property taxes on that parcel.
Thank you for your assistance.
Steven J. Pierce
Center for Statistical Training & Consulting (CSTAT)
Michigan State University
Web: http://www.cstat.msu.edu
-----Original Message-----
From: Roger Bivand [mailto:Roger.Bivand at nhh.no]
Sent: Monday, June 22, 2009 2:59 AM
To: Steven J. Pierce
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Problem reading in polygon shapefile with
readShapePoly
On Mon, 22 Jun 2009, Steven J. Pierce wrote:
> Hi folks,
>
> I've been trying to read a polygon shapefile I received from an ArcGIS
> user into R with the tools provided in the maptools and rgdal
> packages. When I tried reading in a file with boundaries for parcels
> of land, I get an error message that reads "Error in from[j]:to[j] :
> NA/NaN argument" and the parcels object is not created. This is
> similar to a situation recently discussed on this list (on June 10th),
> so I tried adding the option delete_null_obj=TRUE to the
> readShapePoly() command, but then just get a different error message
> after doing that and the object still does not get created. Based on
> the mention of rings not being closed in the new error message, I took
> a guess that maybe adding force_ring = TRUE would help, but that also
> has not solved the problem. The parcels object is still not being created.
>
> The shapefile seems to work just fine in ArcGIS and as part of a
> published map I can open in ArcReader. In fact, I've got access to
> couple different parcel shapefiles and all of them generate this
> error, whereas all the other polygon shapefiles I've got access to
> (e.g., Census block boundaries) work just fine. I'd really like to get
> this parcel shapefile imported into R and stored as a
> SpatialPolygonsDataFrame object, so I would appreciate any insight you
> can provide on what's going wrong or suggestions on how to diagnose
> and fix the problem, whether it that requires changing my code, or even
modifying the shapefile back in ArcGIS.
>
> Below is an excerpt from my console, with traceback() and
> sessionInfo() output at the end.
Thanks for a full report. What this is saying is that there are both NULL
geometries that are being removed, and in addition there appear to be
non-closed rings even when the force_ring policy is being used.
Without access to the file, I'm afraid that finding a solution is not
feasible. If you can zip the parcel_mast_20050802.* files (and possibly
others with the same problem) and put them online, and send me the URL
off-list, or attach the zip archive if < 2MB, I'll have a look. I have an
as-yet unreleased version of maptools that may be able to handle this case,
but seeing the offending files would be very useful.
Hope this helps,
Roger
>
>> parcels <- readShapePoly("Parcels
>> Shapefile/parcel_mast_20050802.shp",
> + proj4string = CRS(parcels.proj4), verbose=TRUE,
> + delete_null_obj=TRUE, force_ring=TRUE)
> Shapefile type: Polygon, (5), # of Shapes: 21842 Error in
> validityMethod(object) : ring not closed In addition: Warning message:
> In .Map2PolyDF(Map, IDs = IDvar, proj4string = proj4string, force_ring
> = force_ring, :
> Null objects with the following indices deleted: 8519, 9096, 9528,
> 13665, 15880
>> traceback()
> 13: stop("ring not closed")
> 12: validityMethod(object)
> 11: identical(x, TRUE)
> 10: anyStrings(validityMethod(object))
> 9: validObject(.Object)
> 8: initialize(value, ...)
> 7: initialize(value, ...)
> 6: new("Polygon", sl, labpt = cG$cents, area = cG$area, hole =
> as.logical(hole),
> ringDir = as.integer(rD))
> 5: Polygon(coords = crds)
> 4: .shp2srsI(shapes[[belongs[[i]][j]]],
> .nParts.shpI(shapes[[belongs[[i]][j]]]),
> force_ring = force_ring)
> 3: .asSpatialPolygonsShapes(Map$Shapes, IDs, proj4string = proj4string,
> force_ring = force_ring)
> 2: .Map2PolyDF(Map, IDs = IDvar, proj4string = proj4string, force_ring
> = force_ring,
> delete_null_obj = delete_null_obj, retrieve_ABS_null =
> retrieve_ABS_null)
> 1: readShapePoly("Parcels Shapefile/parcel_mast_20050802.shp",
> proj4string = CRS(parcels.proj4),
> verbose = TRUE, delete_null_obj = TRUE, force_ring = TRUE)
>> is.projected(parcels)
> Error in extends(class(sd), "Spatial") : object 'parcels' not found
>>
>> sessionInfo()
> R version 2.9.0 (2009-04-17)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
> States.1252;LC_MONETARY=English_United
> States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] spatstat_1.15-3 deldir_0.0-8 gpclib_1.4-4 mgcv_1.5-5
> [5] rgdal_0.6-8 maptools_0.7-23 sp_0.9-37 foreign_0.8-36
>
> loaded via a namespace (and not attached):
> [1] grid_2.9.0 lattice_0.17-25 nlme_3.1-92
>
>
> Steven J. Pierce
> Center for Statistical Training & Consulting (CSTAT) Michigan State
> University
> Web: http://www.cstat.msu.edu
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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