[R-sig-Geo] random spatial distribution code / shapefiles with more than one polygon

[Abby Rudolph]

Hi Marcelino,

I think what I would like to do is convert the shapefile comprised of multiple polygons into a single object of class owin.
However, when I use the code as(NYCBorough, "owin")
I get the following error:
Error in owin(poly = opls): Polygon data contain overlaps between polygons.

Is there a way to get around this?

Regarding the miles part, I need to convert both to miles for a calculation that comes later on in my code.

Thanks,
Abby

-----Original Message-----
From: Marcelino de la Cruz [mailto:marcelino.delacruz at upm.es] 
Sent: Thursday, March 07, 2013 2:53 AM
To: Abby Rudolph
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] random spatial distribution code / shapefiles with more than one polygon

Hi,
My first advise is to point you to the vignette about "handling shape les in the spatstat package" .Go to section 3.2.5 "Objects of class SpatialPolygons".

In short, you could do something like this:

BaltCity2<-readShapePoly("baltcity_mdsp83m.shp")
play2<-ppp(play$X_m,play$Y_m, window=as.owin(BaltCity2))

...if you want to combine all polygons inside baltcity_mdsp83m.shp in just one window. See the section section 3.2.5 "Objects of class SpatialPolygons" for the other option (having different windows from each indiviual polygon).

As you can see, I also removed the "meter to miles" conversion because it's unnecesary  if both the shapes and the  play$X_m and play$Y_m coordinates are in the same units (metres)


HTH,
Marcelino


El 07/03/2013 4:51, Abby Rudolph escribió:
> Hi,
>
> I was given the following code to determine whether the spatial 
> distribution of attribute=1 is significantly different from the spatial distribution of attribute=0 given that the points (overall) are not randomly distributed in space.
>
>
> BaltCity2<-readShapePoly("baltcity_mdsp83m.shp")
> summary(BaltCity2)
>
> x<-BaltCity2 at polygons[[1]]@Polygons
> x
> summary(x)
> # this shows that the baltimore polygon is a single shape # the 
> Baltimore shapefile is projected in meters, so the following code 
> converts meters to miles
> coords<-slot(x[[1]],"coords")
> coords<-coords[167:1,]
> coords<-coords[-1,]
> coords.mi<-coords*0.000621371192
> summary(coords.mi)
>
> play2<-ppp(play$X_m*0.000621371192,play$Y_m*0.000621371192,
>             window=owin(poly=list(x=coords.mi[,1],y=coords.mi[,2])))
>
> I am able to get the program to work when I use the Baltimore shapefile because the polygon is a single shape.
>
> However, when I apply this code to data I have collected from New York City, I have problems.  I am wondering if the problem is the result of the NYC shapefile having multiple polygons and water gaps.  I have been attempting this code with two different NYC shapefiles.  The first is from NYC Bytes.  I brought it into ArcMap and projected it into meters.
> When I plot the following in R, the points from my data file map onto 
> the shapefile and both are in meters
>
> plot(NYCBB)
> u<-RDS_TSO$RDS_new==0
> points(RDS_TSOm$X_m[u],RDS_TSOm$Y_m[u],pch=5,col="yellow")
> points(RDS_TSOm$X_m[!u],RDS_TSOm$Y_m[!u],pch=5,col="red")
>
> However, when I use the code below to convert the meters to miles:
>
> x<-NYBorough at polygons[[1]]@Polygons
> x
> # there are 8910 rows in x
> summary(x)
> # x is 5 different shapes compared to one like in the baltimore example above
>       Length Class   Mode
> [1,] 1      Polygon S4
> [2,] 1      Polygon S4
> [3,] 1      Polygon S4
> [4,] 1      Polygon S4
> [5,] 1      Polygon S4
>
> coords<-slot(x[[1]],"coords")
> coords
> # Now there are only 47 rows - with the baltimore data there was no 
> reduction - I suspect it is being limited in some way 
> coords<-coords[19:1,] coords coords<-coords[-1,] coords
> NYBorough.coords.mi<-coords*0.000621371192
> NYBorough.coords.mi
>
> then, when I try to convert my data x,y coordinates from meters to 
> miles and create a shape with the NYBorough owin, I get a warning that 
> all of the points fall outside the window
>
> RDS<-ppp(RDS_TSOm$X_m[u]*0.000621371192,RDS_TSOm$Y_m[u]*0.000621371192,
>             
> window=owin(poly=list(x=NYBorough.coords.mi[,1],y=NYBorough.coords.mi[
> ,2])))
>
> I received the following warning:
> Warning message:
> In ppp(RDS_TSOm$X_m[u] * 0.000621371192, RDS_TSOm$Y_m[u] * 0.000621371192,  :
>    217 points were rejected as lying outside the specified window
>
> This is puzzling, because when they are plotted using the following 
> code, they overlap
>
> plot(NYBorough)
> points(RDS_TSOm$X_m[u],RDS_TSOm$Y_m[u],pch=5,col="yellow")
>
> and there are only 217 points in this group, so it appears that they are all rejected.
>
> I've also tried downloading the TIGER files from the US Census and merging all 5 boroughs, projecting the data in arcmap into meters and opening in R.
> I get the same warning with both shapefiles.
>
> Does it seem like the error is resulting from the shapefile I am using or some other error in the code?
>
> Any help would be very much appreciated.
>
> Thanks!
>
> Abby
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


--
Marcelino de la Cruz Rot
Depto. de Biología Vegetal
Universidad Politécnica de Madrid
Madrid, España