7

I am trying to break out a single feature of a Spatialpolygon object. That is e.g. plot a single polygon from a shapefile that contains multiple. Any suggestions how to do this?

>#Import shapefile, generate spatialpolygondataframe
>vectormap<-readOGR(dsn="QGIS export", layer="section4") 
>
>OGR data source with driver: ESRI Shapefile 
>Source: "QGIS export", layer: "section4"
>with 4 features
>It has 6 fields
>
>
># Change vectormap CRS to Equirectangular Carree projection
>vectormap2 <- spTransform(vectormap, CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 >+x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

> #Summary of vectormap2
>Object of class SpatialPolygonsDataFrame
>Coordinates:
>       min      max
>x  1467755  5288689
>y -6018287 -3018986
>Is projected: TRUE 
>proj4string :........
>Data attributes:
>  Layer                       SubClasses ExtendedEn Linetype EntityHand Text  
> Layer 1:4   AcDbEntity:AcDb2dPolyline:3    NA's:4     NA's:4   41:1    NA's:4   
>             AcDbEntity:AcDbSpline    :1                        5B:1               
>                                                                7B:1               
>                                                                94:1               
> coordinates(vectormap2)
>     [,1]     [,2]
>0 4122892 -3791580
>1 4661889 -4500891
>2 4500413 -4079798
>3 3352295 -4392646

Happy for any input!

Side Note. I can also get the shapefile in a format giving the features in separate layers.

> summary(vectormap)
>Object of class SpatialPolygonsDataFrame
>Coordinates:
>           min      max
>x  -0.01383442 457.4256
>y 125.27349310 484.3477
>Is projected: FALSE 
>proj4string :............
>Data attributes:
>         Layer                           SubClasses ExtendedEn Linetype >EntityHand   Text  
> big        :1   AcDbEntity:AcDbBlockReference:4    NA's:4     NA's:4   42:1       >NA's:4  
> fullsection:1                                                          48:1               
> medium     :1                                                          67:1               
> small      :1 
84:1

However when I do this, I cant change the projection.

vectormap2 <- spTransform(vectormap, CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
non finite transformation detected:
[1] 187.4585 468.4831      Inf      Inf
Error in .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  : 
  failure in Polygons 1 Polygon 1 points 1
In addition: Warning message:
In .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  1808 projected point(s) not finite

1 Answer 1

10

You need to select which feature will be plotted. For this purpose, select the observation as the same way as you select a row or more from a data.frame object:

plot(shp)

enter image description here

plot(shp[1,])

enter image description here

plot(shp[2:3,])

enter image description here

plot(shp[shp$Nombre == 'Estanque',])

enter image description here


This also work to extract features to a new polygon:

shp2 <- shp[1,]

shp2

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 317185.2, 317341.5, 6601234, 6601492  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 3
names       : id,    Nombre,   Sup 
min values  :  1, Cuartel 1, 1.937 
max values  :  1, Cuartel 1, 1.937 
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.