6

I am to map a plot of an area and then plot points over it (given by Easting and Northing). However, I found out that there are certain mistakes in the data entered for Easting and Northing resulting in some outlying points, and ggplot (which I am using) tends to reshape the map in order to fit the outlying points.

Is there any quick way to fix this issue? I do not know of any except to identify and correct the points manually. I am using shapefiles from the Ordnance survey.

Here is the screen-shot. The image does not appear to be skewed. However, upon comparison with other images, the differences become apparent.

**

EDIT

**

The code I'm using is as follows.

countyRegion<- readShapePoly(file.choose())
norfolkCounty<- countyRegion[countyRegion$NAME=="Norfolk County",]
#convert shape file into data that can be plotted on graph
gpclibPermit()
norfolk<- fortify(norfolkCounty,region="NAME")
Norfolk<- merge(norfolk, norfolkCounty@data,by.x="id",by.y="NAME")

#eit data
jul2012<- read.table(file.choose(),header=TRUE,as.is=TRUE,blank.lines.skip=FALSE,sep=",")
jul2012$Crime.type<- factor(jul2012$Crime.type)
jul2012<- jul2012[jul2012$Crime.type!="",]
jul2012$Crime.type<- factor(jul2012$Crime.type)
levels(jul2012$Crime.type)<- c("Anti-social behaviour","Theft/burglary","Criminal damage/arson","Drugs","other crime","Theft/burglary","Public disorder and weapons","Theft/burglary","Vehicle crime","Violence and Sexual Offences")

#plot
jul12Map<- ggplot(data=Norfolk,aes(long,lat)) + geom_polygon() + geom_point(data=jul2012,aes(Easting,Northing,colour=as.factor(Crime.type)),alpha=0.6) + scale_colour_brewer(palette="Set1",name="Category") + theme_bw() + scale_x_continuous(breaks=0,labels="") + scale_y_continuous(breaks=0,labels="") + theme(axis.ticks=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank())

The shp file and the crime stats file are linked, as well as the dbf, shx, and prj files. enter image description here enter image description here

  • 1
    Agreed. It's impossible to tell whether it's skewed and by how much from this screenshot. I know the OS files are large, but if you can supply a link and a section of your data that contains the outlier points we may be able to help. – SlowLearner Aug 26 '13 at 20:33
  • 1
    Without seeing the code that produced this, we can't help. Did you have +coord_equal() in your ggplot to make sure the aspect ratio is unity? Is that what you mean by 'skewed'? Because that's 'squashed' rather than 'skewed'. Removing the outlying points is a separate question. – Spacedman Aug 26 '13 at 21:29
  • You're right, I should have been clearer. I have now edited the question to provide images, the code I use, and links to the shape and data files. Apologies for the ambiguity. – info_seekeR Aug 27 '13 at 3:10
  • Much better, thanks for that. Unfortunately as well as the county_region.shp file we need the .avl, .dbf, .prj and .shx files as well. It is the combination of these that makes a useable shapefile. – SlowLearner Aug 27 '13 at 6:31
  • Oh, pardon me again. I recently started working with GIS in spare time, so I am relatively new to the field. I have uploaded .dbf, .prj, .shx files, but I hadn't received any .avl files in my download from OS... – info_seekeR Aug 27 '13 at 7:49
5

I'm assuming that you are using the Meridian 2 data from the Ordnance Survey and have used that and your data to create the chart below. I have stripped away the parts of your code that seemed to me to be redundant. Basically as far as I can see you want to automatically remove points that do not fall within the area in question, in this case the county of Norfolk.

The key issue (it seems to me) is to create a SpatialPoints object and use the over function from the sp package to check whether the points lie within the polygon of interest. This post provides useful guidance. If the points lie outside the polygon, you end up with NA values, which can be found easily enough. (EDIT Incidentally, the points that lie outside the Norfolk polygon do refer to valid addresses - they are not simple transposition mistakes of Eastings or Northings.)

Here's what the map looks like with outliers:

norfolk bad

Here's what the map looks like with outliers removed:

norfolk good

Here is the code used to create the maps. I haven't had much experience with over and similar functions so you probably want to experiment to ensure the code is sound. I checked the 'bad' points and they seemed in most cases to correspond with the outliers on the map.

require(rgdal)
require(ggplot2)

work.dir <- "your_dir_here_no_trailing_slash"

countyRegion <- readOGR(work.dir, layer = "county_region")
norfolkCounty<- countyRegion[countyRegion$NAME == "NORFOLK COUNTY",]

norfolk <- fortify(norfolkCounty, region = "NAME")

jul2012<- read.csv(paste0(work.dir, "/2012-07-norfolk-street.csv"),
                     header = TRUE)
names(jul2012)[which(names(jul2012) == 'Easting')] <- 'long'
names(jul2012)[which(names(jul2012) == 'Northing')] <- 'lat'

ggplot(data = norfolk, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "black", fill = "white") +
    scale_x_continuous(breaks = seq(500000, 700000, 5000)) + # makes it easier to
    scale_y_continuous(breaks = seq(240000, 350000, 5000)) + # locate points on map
    coord_equal() +
    theme_minimal() +
    geom_point(data = jul2012,
               size = 2,
               pch = 21,
               alpha = 0.6,
               colour = "white",
               aes(x = long, y = lat, group = NULL,
               fill = as.factor(Crime.type))
               ) +
    theme(legend.position = "none")

# now remove points that are not contained within polygon of interest
pts <- SpatialPoints(jul2012[, c('long','lat')], proj4string = norfolkCounty@proj4string)
pts.over <- over(pts, norfolkCounty)
print(str(pts.over)) # check
pts.bad <- rownames(pts.over[is.na(pts.over$IDENTIFI0), ])
print(jul2012[pts.bad, ]) # how many do we have
jul2012.good <- jul2012[!(rownames(jul2012) %in% pts.bad), ]

ggplot(data = norfolk, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "black", fill = "white") +
    scale_x_continuous(breaks = seq(500000, 700000, 5000)) +
    scale_y_continuous(breaks = seq(240000, 350000, 5000)) +
    coord_equal() +
    theme_minimal() +
    geom_point(data = jul2012.good,
               size = 2,
               pch = 21,
               alpha = 0.6,
               colour = "white",
               aes(x = long, y = lat, group = NULL,
               fill = as.factor(Crime.type))
               ) +
    theme(legend.position = "none")
  • Thanks so much for kindly experimenting with the code and presenting a solution. I take your point that you do not have a lot of experience with functions such as over. I will experiment with it, and post my results here too. I'm not sure that I will be able to add my answer if I accept yours first, so I will wait a while, post my results, and then accept your answer. Thanks for the link too! – info_seekeR Aug 28 '13 at 10:19

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.