Skip to content

Intersect & Plot a Polygon and wkbLineString

5 messages · Jim Burke, Roger Bivand

Original
Jim Burke
# 
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 **

b)How to plot the result 

c)How to display the street name on the plot (that's FNAME from the 
major roads.dbf file).
*


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
Roger Bivand
#
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).
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).
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?
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

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jim Burke
    

    

      
      #
    

  


  

      
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.

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.

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

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

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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Tue, 17 Feb 2009, Jim Burke wrote:

        

            

      
      
      
    

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

            

      
      
      
    

        
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.

            

      
      
      
    

        
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?

            

      
      
      
    

        
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