Tell me more ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 2 down vote favorite

I am using the PBSmapping package in R to model data from diva-gis. I can read in the shapefiles and plot them with no issue.

library(PBSmapping)
df = importShapefile("/path/df.shp")
plotPolys(df)

But the problem is I can't figure out how to subset the PolySet data. For example, for Bhutan's administrative areas file, runnning str(df) shows:

Classes ‘PolySet’ and 'data.frame': 59294 obs. of  5 variables:
 $ PID: int  1 1 1 1 1 1 1 1 1 1 ...
 $ SID: int  1 1 1 1 1 1 1 1 1 1 ...
 $ POS: int  1 2 3 4 5 6 7 8 9 10 ...
 $ X  : num  91 91 91 91 91 ...
 $ Y  : num  27.2 27.2 27.2 27.2 27.2 ...
 - attr(*, "PolyData")=Classes ‘PolyData’ and 'data.frame': 201 obs. of  19 variables:
  ..$ PID       : int  1 2 3 4 5 6 7 8 9 10 ...
  ..$ ID_0      : int  35 35 35 35 35 35 35 35 35 35 ...
  ..$ ISO       : Factor w/ 1 level "BTN": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ NAME_0    : Factor w/ 1 level "Bhutan": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ ID_1      : int  486 486 486 486 486 486 486 486 472 472 ...
  ..$ NAME_1    : Factor w/ 20 levels "Bumthang","Chhukha",..: 15 15 15 15 15 15 15 15 1 1 ...
  ..$ ID_2      : int  9280 9281 9282 9283 9284 9285 9286 9287 9141 9142 ...
  ..$ NAME_2    : Factor w/ 201 levels "Aemeri","Athang",..: 8 19 73 126 129 138 163 186 4 29 ...
  ..$ VARNAME_2 : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
  ..$ NL_NAME_2 : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
  ..$ HASC_2    : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
  ..$ CC_2      : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
  ..$ TYPE_2    : Factor w/ 1 level "Unknown": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ ENGTYPE_2 : Factor w/ 1 level "Unknown": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ VALIDFR_2 : Factor w/ 1 level "Unknown": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ VALIDTO_2 : Factor w/ 1 level "Unknown": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ REMARKS_2 : Factor w/ 1 level "Second level is probably 3rd level": 1 1 1 1 1 1 1 1 1 1 ...
  ..$ Shape_Leng: num  0.82 0.706 0.433 1.053 0.751 ...
  ..$ Shape_Area: num  0.01852 0.01736 0.00877 0.04298 0.01871 ...
 - attr(*, "parent.child")= num  1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "shpType")= int 5
 - attr(*, "prj")= chr "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"D"| __truncated__
 - attr(*, "projection")= chr "LL"

How can I subset out PolySet data where "Name_1" is "Bumthang"? Incidentally, in maptools, this can be done by:

library(maptools)
df = readShapeSpatial("/path/df.shp")
df.data = as.data.frame(df) #For understanding the data structure
df.Bumthang = admin.mpt[as.character(admin.mpt$NAME_1)=='Bumthang',]
share|improve this question
I can't get the code to run... can you post some straightforward instructions on how to get "Bhutan's administrative area" so I can duplicate your "df" ? – J. Winchester Feb 22 '12 at 8:06
Sure. Go to: diva-gis.org/gdata, and select "Bhutan" (or any other country) and "Administrative Areas" and hit "OK". The next page will show a map of Bhutan, with a red "Download" link at the bottom. This will download the "BTN_adm.zip" file. You can use the "BTN_adm2.shp" file in there. – user3671 Feb 22 '12 at 11:32
Does it have have to be a Polyset? Won't it be easier to subset it first via some other module, write to file, and then plot the resulting subset via PBSmapping? – R.K. Mar 1 '12 at 5:51

1 Answer

up vote 2 down vote

What about using a using another library to subset the file? You could subset the data using rgdal, save to file, import the shapefile containing the subset and plot it using PBSMapping.

library(rgdal)
library(PBSmapping)

df <-  readOGR(".","df")
subset <- df[df$NAME_1=="Bumthang",]
writeOGR(subset, ".", "bumthang", driver="ESRI Shapefile")
bum = importShapefile("bumthang.shp")
plotPolys(bum, projection=TRUE)

enter image description here

Or you could do away with PBSmapping altogether.

 plot(subset, axes=TRUE)

enter image description here

share|improve this answer
I was trying to find something within PBSmapping. One alternative is to use PolySet2SpatialPolygons in the maptools library, and then subset it. Youre rgdal pointer is nice. Many thanks, R.K. – user3671 Mar 4 '12 at 0:38
It's just that using the PolySet can be slow. I've tried it on my machine for comparison and computer froze for a moment when I tried the PBSmapping solution. – R.K. Mar 4 '12 at 2:43

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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