Выбор правильного значения для proj4string для чтения шейп-файла в R?

У меня есть шейп-файл полигонов и другой файл CSV, который содержит список точек в виде (Lat, Lng) пар.

Я хочу проверить для каждой пары (lat, lng) из файла CSV, в какой полигон он попадает внутрь ..

Шейп-файл проецируется, и файл proj выглядит так:

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

Мой план следующий:

1. Прочтите шейп- файл, используя readShapePolyфункцию из MapToolsпакета R.
2. Считать координаты точек из CSV-файла в информационный кадр и преобразовать его в SpatialPointsDataFrame
3. Используйте overфункцию, чтобы определить, в какой полигон он попадает.

Для этого мне нужно указать proj4stringвремя загрузки шейп-файла на шаге 1, а также преобразовать координаты из файла CSV в ту же систему проекции, используя spTransformфункцию, прежде чем применять overфункцию на шаге 3, так как для этого необходимо, чтобы точки и многоугольники быть под той же системой проекции.

Любая идея о том, что должно правильное значение для содержимого файла proj, показанного выше?

Proj4string является допустимой строкой PROJ4 crs.

см. Как я могу получить строку proj4 или код EPSG из файла shapefile .prj? и Shapefile PRJ для таблицы поиска PostGIS SRID?

короче говоря:

• Вы можете использовать gdalinfo как в первой ссылке или привязки GDAL Python как во второй ссылке.

Или

• перейти к Prj2EPSG (простой сервис для преобразования хорошо известной информации о проекции текста из файлов .prj в стандартные коды EPSG)
• Введите содержание вашего файла prj

• результат EPSG: 27700, поэтому первая версия строки PROJ4

  " + init = epsg: 27700 "

`Или

• Перейдите в Spatial Reference и введите код -> результат EPSG Projection 27700 - OSGB 1936 / British National Grid

• нажмите на Proj4, и полная строка PROJ4:

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

Вот очень удобный сайт для получения кода EPSG для данного прогноза. В вашем случае прогноз "EPSG: 27700". Если у вас есть проекции, определенные для шейп-файла, вы можете назначить проекцию при создании SpatialPointsDataFrame, а затем использовать определение проекции из импортированного шейп-файла. Использование «readOGR» из пакета rgdal сохранит определения проекции. Вот пример назначения и извлечения координатных строк в данных класса sp.

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