Hoe selecteer ik een subset van PolySet-gegevens in R?

Ik gebruik het PBSmapping -pakket in R om gegevens te modelleren van diva-gis . Ik kan in de shapefiles lezen en ze plotten zonder problemen.

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

Maar het probleem is dat ik niet weet hoe ik de PolySet-gegevens moet subsetten. In het bestand met administratieve zones van Bhutan wordt bijvoorbeeld str (df) weergegeven:

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"

Hoe kan ik PolySet-gegevens uitschakelen waar 'Name_1' 'Bumthang' is? Overigens kan dit in maptools worden gedaan door:

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',]
4
Moet het een Polyset zijn? Zal het niet eenvoudiger zijn om het eerst via een andere module onder te verdelen, naar bestand te schrijven en vervolgens de resulterende subset uit te zetten via PBSmapping?
toegevoegd de auteur shsteimer, de bron
Ik kan de code niet laten uitvoeren ... kun je enkele duidelijke instructies plaatsen over hoe je "het administratieve gedeelte van Bhutan" kunt krijgen, zodat ik je "df" kan dupliceren?
toegevoegd de auteur Davide Vosti, de bron
Zeker. Ga naar: diva-gis.org/gdata en selecteer "Bhutan" (of een ander land) en "Administratieve gebieden" en druk op "OK". De volgende pagina toont een kaart van Bhutan, met onderaan een rode "Download" -link. Hiermee wordt het bestand "BTN_adm.zip" gedownload. U kunt het bestand "BTN_adm2.shp" daar gebruiken.
toegevoegd de auteur Chris Canal, de bron

1 antwoord

Hoe zit het met het gebruik van een andere bibliotheek om het bestand subset te maken? U kunt de gegevens subset met rgdal, opslaan in bestand, de shapefile importeren die de subset bevat en deze plotten met 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

Of je zou PBSmapping helemaal kunnen schrappen.

 plot(subset, axes=TRUE)

enter image description here

4
toegevoegd
Het is alleen dat het gebruik van de PolySet langzaam kan zijn. Ik heb het op mijn computer geprobeerd voor vergelijking en de computer bevroor voor een moment toen ik de PBSmapping-oplossing probeerde.
toegevoegd de auteur shsteimer, de bron
Ik probeerde iets te vinden binnen PBSmapping. Een alternatief is om PolySet2SpatialPolygons in de maptools-bibliotheek te gebruiken en deze vervolgens te subsetten. Je rgdal pointer is leuk. Hartelijk dank, R.K.
toegevoegd de auteur Chris Canal, de bron