Building elements from shapefiles in R

Question

stackoverflow community

I have two forms from the municipalities of Japan. I use R to create separate plots for each of the municipalities. I can do this work with one shapefile, but identical syntax fails with another. Here is the code and URLs for the data:

library(sp) library(maptools) # Map A - this one works # Please download: http://www.filefactory.com/file/z26nirxoz53/n/JPN_adm_zip # Enter your path for readShapePoly japanMapA = readShapePoly("JPN_adm/JPN_adm2") names(japanMapA) plot(japanMapA[japanMapA$ID_2 == 1199,]) # Map B - this one doesn't work # Please download: http://geocommons.com/overlays/173340.zip # Again, enter your path for readShapePoly japanMapB = readShapePoly("japan_ver71") names(japanMapB) plot(japanMapB[japanMapB$JCODE == 45382,]) 

The error it throws is:

 Error in plot(japanMapB[japanMapB$JCODE == 45382, ]) : error in evaluating the argument 'x' in selecting a method for function 'plot': Error in japanMapB[japanMapB$JCODE == 45382, ] : NAs not permitted in row index 

I do not know how to do this in order to remove NA in this case, so I cannot build individual elements.

I would thank you for your help: I ​​banged my head on the wall for a few seconds!

Answer 1

I think it comes down to a problem with the appropriate setting.

Your japanMapB object consists of several metadata and a series of polygons for each form stored in japanMapB@polygons . So you have:

 > length(japanMapB$JCODE) #[1] 1902 > length( japanMapB@polygons ) #[1] 1902 

As @PaulHiemstra notes, you have multiple NA values ​​in your JCODE variable

 > table(is.na(japanMapB$JCODE)) #FALSE TRUE # 1894 8 

This means that you get NA results when you try to index the municipalities that you want to build.

 > table(japanMapB$JCODE==45382,useNA="always") #FALSE TRUE <NA> # 1893 1 8 

Wrapping in which solves the following:

 which(japanMapB$JCODE == 45382) #[1] 1802 #You will now select to plot only the 1802th polygon stored in the data object plot(japanMapB[which(japanMapB$JCODE == 45382),]) 

Answer 2

What can happen is that there are municipalities / polygons with the identifier NA . For instance:

 id = c(1,2,3,NA) id == 2 # [1] FALSE TRUE FALSE NA 

Please note that in this comparison, the NA values ​​in id do not lead to FALSE , but to NA . Using id == 2 for a subset leads to the problem you see. @thelatemail ends the call in which , which translates the logical vector into a set of indices, where the logical vector is TRUE :

 which(id == 2) # [1] 3 

Note that calling which excludes NA . This index vector can be safely used for a subset. In general, you can use na.omit to remove NA :

 l[l == 2] # [1] 2 NA l[na.omit(l == 2)] # [1] 2