[R-sig-Geo] Intersect & Plot a Polygon and wkbLineString

[Roger Bivand]

On Tue, 17 Feb 2009, Jim Burke wrote:

> Sorry, we seem to have misunderstood each other.  Basically I have a 
> county-wide set of roads called "major roads".  Then originally I had an 
> identically sized county-wide set of polygons of which I select (via R) 
> a smaller than county-wide subset of adjacent polygons and plot that 
> subset.
>
> The "add = TRUE" in plot operations suggestion sounds interesting. Your 
> last paragraph is indeed all I am trying to achieve.  I need a GIS 
> program to do this? I thought that R might be able to clip the major 
> roads file (via an intersect) to conform to the polygon shapes.

Thanks, that is clearer. There is some support for polygon clipping in the 
gpclib package, that other packages also use, but not for lines.

>
> I don't know, maybe I am approaching this the wrong way (i.e. trying to 
> find the SQL style intersection of these two blobs). Perhaps another 
> approach may be better. I think I will do the following.
>
> 1. Plot the major roads and try to get the street names on them. No 
> polygon here.

OK, perhaps using xlim= and ylim= to "zoom" in on the view, and locator() 
to choose the label points. Maybe with locator(2) and a bit of 
trigonometry, you could get a label point and an angle. But 30000 is a 
lot for heads-up digitizing.

You might start with some ideas like:

library(maptools)
fylk <- readShapeSpatial(system.file("shapes/fylk-val.shp", package =
   "maptools"))
library(spatstat)
y <- as(fylk, "SpatialLines")
z <- as(y, "psp")
plot(z)
zz <- midpoints.psp(z)
plot(zz, add=TRUE, col="green", pch=3, cex=0.2)
df <- data.frame(x=zz$x, y=zz$y, seg=z$mark)
table(df$seg)

gives a range of segment midpoints - where the discrete values of df$seg 
would indicate which Lines object (FNAME) they correspond to. If there is 
only one, you're OK, if more, maybe choose a central one. They you could 
use text() with the coordinates of the midpoints, one for each FNAME 
object, which would show you where you are roughly. Choosing the midpoints 
ought to keep you away from the junctions. If the roads are gridded, the 
angles will be regular, otherwise that'll need some attention.

>
> 2. Overlay what's going to be the large major roads over the smaller 
> polygon set. The "add = TRUE" suggestion. Thanks.
>

That means working out which of your 30000 constitute the large major 
roads - if you have an attribute, you can perhaps subset to those roads 
first?

> 3. Figure out how to trim that major roads set.
>

You could coerce the polygons to SpatialLines, and go out to psp in 
spatstat as above. Then you have crossing.psp which shows where two psp 
objects cross, but as far as I can see no indication of which segments 
cross, so you'd need to subset those.

zzz <- psp(20000, 6600000, 1060000, 7700000, window=zz$window)
plot(zzz, add=TRUE, col="red")
crossing.psp(zzz, z)
plot(crossing.psp(zzz, z), add=TRUE, col="blue", pch=3)

Use unionSpatialPolygons() in maptools to make a cookie cutter from the 
external boundary of the polygons, then you only need to loop over the 
roads and find out which segment needs the point inserted, cutting it 
there. Another possibility is to fake it by taking the complement of the 
outer edge of the polygons, and over-plot that, painting over the roads.

Hope this helps,

Roger


> I have your book "Applied Spatial Data Analysis with R (Use R)" on order 
> and it should arrive later today. Then also I am looking forward to a 
> good read of "Lattice Multivariate Visualization".
>
> Thanks,
> Jim Burke
>
>
>
> Roger Bivand wrote:
>> On Mon, 16 Feb 2009, Jim Burke wrote:
>> 
>> (Just stick to plain text, but don't use fonts or bold or stuff - keep it 
>> simple).
>> 
>>> I have a nice solid SpatialPologonsDataFrame and a new line file 
>>> (wkbLineString) called "major roads.shp".
>>> 
>>> MY QUESTIONS.
>>> 
>>> *I would appreciate examples of*
>>> 
>>> *a)How to merge the intersection of a **SpatialPologonsDataFrame and a 
>>> wkbLineString **
>>> 
>> 
>> Why, and what do you mean? You have some 30000 road lines, what do you want 
>> to know? How much of which road is in which polygon? Again, you need to 
>> think through your workflow. My inclination would be to rasterise the lines 
>> in a GIS, and overlay the raster cell centre points on the polygons, but 
>> maybe you want line length, or to retain the line IDs (all 30000 of them).
>> 
>> 
>>> b)How to plot the result
>> 
>> You can of course overplot lines on polygons, with add=TRUE in plot() 
>> methods and sp.layout= in spplot() methods. But plotting the result of the 
>> "intersection" means what?
>> 
>>> c)How to display the street name on the plot (that's FNAME from the major 
>>> roads.dbf file).
>> 
>> 30000 names? Do you have a very large screen? Lines objects have IDs, but 
>> no label point, so there is no easy way of doing it, even label alignment 
>> is hard.
>> 
>> I can't actually see the statistical question here, could you please make 
>> it clearer?
>> 
>> Do you simply want to plot a subset of the road lines on your solid 
>> polygons, adding names? Then you need to find out how to subset them, 
>> perhaps by cookie-cutting using the union of your polygons, add label 
>> points and rotation angles (probably by hand), and off you go, but these 
>> are essentially GIS operations.
>> 
>> Roger
>> 
>>> *
>>> 
>>> 
>>> ADDITIONAL INFORMATION.
>>> 
>>> The SpatialPologonsDataFrame is a garden variety collection of 
>>> geographical polygons.
>>>> tx1_sp <- readShapePoly("precinct08.shp", IDvar="PCT", 
>>> proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))
>>> 
>>> 
>>> 
>>> The new file to intersect merge is major roads line map that should 
>>> overlay part of my SpatialPolygons. It does not behave like a Shapfile.
>>> 
>>>> main_roads <-readOGR("main roads.shp", layer = "main roads")
>>> OGR data source with driver: ESRI Shapefile
>>> Source: "main roads.shp", layer: "main roads"
>>> with  29161  rows and  33  columns
>>> Feature type: wkbLineString with 2 dimensions
>>> 
>>> 
>>> Thanks,
>>> Jim Burke
>>> 
>>> _______________________________________________
>>> 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