spatiallyConfused spatiallyConfused - 4 months ago 28
R Question

Plotting elements from shapefiles in R

stackoverflow community,

I have two shapefiles of the municipalities of Japan. I'm using R to create separate plots for each of the municipalities. I can make this work with one shapefile, but the identical syntax fails with the other. Here are the code and the URLs for the data:


# Map A - this one works
# Please download:

# Enter your path for readShapePoly
japanMapA = readShapePoly("JPN_adm/JPN_adm2")

plot(japanMapA[japanMapA$ID_2 == 1199,])

# Map B - this one doesn't work
# Please download:

# Again, enter your path for readShapePoly
japanMapB = readShapePoly("japan_ver71")

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 don't know how to go about removing the NAs in this case, so I am unable to plot the individual elements.

Would greatly appreciate your help: have been banging my head against the wall for a while on this one!


I think this comes down to an issue with appropriate sub-setting.

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

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

As @PaulHiemstra notes though, you have some NA values in your JCODE variable

> table($JCODE))

# 1894     8 

Which means you get NA results when trying to index the municipalities you want to plot.

> table(japanMapB$JCODE==45382,useNA="always")

# 1893     1     8 

Wrapping in which solves this:

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),])