Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am having a shape file 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 shape file 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 fllows:

  1. Read the shape file 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 shape file 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 ?

share|improve this question
    
If your shapefile(s) already have the projection defined use "readOGR" in the rgdal package. This package is a wrapper for GDAL and really supersedes the shapefile read/write functionality in maptools. This function handles all topology types and retains the projection information. –  Jeffrey Evans Jun 27 '13 at 16:52
    
When I try loadign the shape file using readOGR function I a always get Cannot open file error –  Moustafa Alzantot Jun 27 '13 at 17:00
    
OK, Now I have been able to read the file using readOGR.. using summary function for the SpatialPolygonDataFrame object gave me the correct value for the proj4string –  Moustafa Alzantot Jun 27 '13 at 17:14
    
Well, without details on how you are using the function it is difficult to help you! Part of the syntax is the directory the data resides in and you do not need the .shp extension in the filename. Something like readOGR(getwd(), "YourShape") should work if you have your work-directory set to the same place your shepfile is. –  Jeffrey Evans Jun 27 '13 at 17:15
    
Thanks @JeffreyEvans, it worked now and I used it to get the proj4string –  Moustafa Alzantot Jun 27 '13 at 17:20
add comment

2 Answers

up vote 4 down vote accepted

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

enter image description here

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

    "+init=epsg:27700"

`Or

enter image description here

  • 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"

share|improve this answer
add comment

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))   
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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