1

I want to plot a figure with my data combined to an shapefile. However, at the moment I am not able to receive output yet since the data from the shapefile covers a larger area than my dataset. My shapefile covers the whole Netherlands (postal codes 1011: 9999) and my dataset only a small part (postal codes 2651: 3344).

In other words, I need to cut my shapefile to ensure both areas have overlap. Doing so I will have no NA's and I can get output. I have read several posts, including this one: How to cut a shape ile?

However, it is still not clear for me how to handle this. Should I use ogr2ogr within R and can I best code this? The explanation on http://www.gdal.org/ogr2ogr.html is quite difficult to follow.


I selected only for the relevant postalcodes, via the comments of Spacedman and Isaac Freitas (both worked). After this I wanted to plot the data using the following code:

plot(studyArea,col=gray(studyArea$data$sumlimitation/max(studyArea$data$sumlimitation)))

However, I still get the same error message as before: Error in gray(studyArea$data$sumBeperking/max(studyArea$data$sumBeperking)) : invalid gray level, must be in [0,1].

What would cause this even though the data is correct now?

  • Do you want to modify the data and then plot this new data or do you want to keep the data unchanged and just plot a subset of this data? Try crop if you want to change the data inside-r.org/packages/cran/raster/docs/crop. – Iris Jun 9 '16 at 13:32
  • Does your shapefile have the postal codes? And you just want the shapefiles with the postalcodes from 2651 to 3344? Then you just subset, like studyArea = netherlands[netherlands$postal >= 2651 & netherlands$postal <= 3344,] - its like taking rows from a dataframe. Here I assume your postal codes are numeric. – Spacedman Jun 9 '16 at 13:52
  • 1
    @user110743 It looks like you've created a second account instead of using your original one...If you use your original account, you'll be able to edit your question without needing it reviewed first and you can comment on answers to your questions, etc. You should also merge your accounts. – Evil Genius Jun 9 '16 at 14:35
  • @Spacedman I now know, what is the last piece of the puzzle. Some numbers between 2651 and 3344 are missing (non existed areas). Therefore, some NA's are still left in the Spatial Polygon Dataframe, which is the reason for the problems with plotting. Is there a way to only use the numbers in my dataset (so skip the non existent ones) or remove the NA's? – Keizer Jun 9 '16 at 15:20
  • Forgot to add: I am not able to remove them with: studyArea$data <- na.omit(studyArea$data) – Keizer Jun 9 '16 at 15:38
1

I answer this assuming you'll be using rgdal to read the shapefile in as a spatial dataframe. Instead of worrying about cutting, you can make a join using the tigris package. You can use the function geo_join() to combine a spatial dataframe and a dataframe. To only include the subset of postal codes from the dataset, set the "how" option to "inner" to get only the data which is in both the shapefile data and the other dataset.

library(rgdal)
library(tigris)

# read in the shapefile
shape_data <- readOGR(dsn="/home/user/file_name.shp", "file_name")
# read in whatever data you have
data <- read.csv("your_data.csv")
# join the two by common id variable
merged_data <- geo_join(spatial_data=shape_data, data_frame=data, by_sp="postal", by_df="postal", how = "inner")
plot(merged_data)

I had considered cutting a shapefile myself before I discovered this strategy. You can then write out the shapefile if you'd like with writeOGR.

EDIT: For your additional question about the gray() function, you need to make sure whatever calculation you are using is between 0 and 1.

  • Did you try gray(na.omit(studyArea$data$sumlimitation)/max(studyArea$data$sumlimitation)). Sorry, can't comment outside my own answer/comments section yet, still new on the site. – Isaac Freitas Jun 9 '16 at 15:59
  • No problem, but yes I tried that, but unfortunately it does not remove the NA's – Keizer Jun 10 '16 at 7:01
  • instead of tigris::geo_join, you could use the more standard merge – Robert Hijmans Jun 12 '16 at 4:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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