Choosing the correct value for proj4string for shapefile reading in R?

Question

I am having a shapefile of polygons and another CSV file which contains a list of points as (Lat, Lng) pairs..

I want to check for each (lat, lng) pair from the CSV file which polygon does it fall inside..

The shapefile is projected and the proj file reads like this:

PROJCS["Transverse_Mercator",GEOGCS["GCS_OSGB 1936",
DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

My plan is as follows:

1. Read the shapefile using the readShapePoly function in the R MapTools package.
2. Read the points coordinates from the CSV file into a dataframe and convert it to SpatialPointsDataFrame
3. Use over function to determine which polygon it falls inside.

In order to do so, I need to specify the proj4string while loading the shapefile in step 1 and also transform the coordinates from the CSV file into the same projection system using spTransform function before applying the over function in step 3 as it requires that the points and polygons must be under the same projection system.

Any idea about what should the correct value for the proj file content shown above ?

Answer 1

The proj4string is a valid PROJ4 crs string.

see How can I get the proj4 string or EPSG code from a shapefile .prj file? and Shapefile PRJ to PostGIS SRID lookup table?

in short:

• You can use gdalinfo as in the first reference or the GDAL Python bindings as in the second reference.

Or

• go to Prj2EPSG (a simple service for converting well-known text projection information from .prj files into standard EPSG codes)
• Enter the content of your prj file

• the result is EPSG:27700 so a first version of the PROJ4 string is

  "+init=epsg:27700"

`Or

• Go to Spatial Reference and enter the code -> result EPSG Projection 27700 - OSGB 1936 / British National Grid

• click on Proj4 and the complete PROJ4 string is:

  "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"

Answer 2

Here is a very handy website for retrieving the EPSG code for a given projection. In your case the projection is "EPSG:27700". If you have projections defined for the shapefile you can assign the projection when you create the SpatialPointsDataFrame and then use the projection definition from your imported shapefile. Using "readOGR" from the rgdal package will retain the projection definitions. Here is an example of assigning and pulling coordinate strings on sp class data.

require(sp)
require(rgdal)

# Use meuse dataset
data(meuse)

# Coerce into SpatialPointsDataframe
coordinates(meuse) <- ~x+y

# Assign projection
proj4string(meuse) <- CRS("+init=epsg:28992")

# Pull some observations and transform to Lat/Long
meuse.geo <- meuse[sample(dim(meuse)[1],10),]
  prj.LatLong <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
    meuse.geo <- spTransform(meuse.geo, prj.LatLong)

# Pull projection string from meuse.geo and use in spTransform
#   to reproject meuse to lat/long  
( prj <- proj4string(meuse.geo) )   
meuse <- spTransform(meuse, CRS(prj))