Escolhendo o valor correto para proj4string para leitura de shapefile em R?

14

Estou tendo um shapefile de polígonos e outro arquivo CSV que contém uma lista de pontos como pares (Lat, Lng) ..

Quero verificar para cada par (lat, lng) do arquivo CSV em que polígono se enquadra ..

O shapefile é projetado e o arquivo proj é assim:

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

Meu plano é o seguinte:

  1. Leia o shapefile usando a readShapePolyfunção no MapToolspacote R.
  2. Leia as coordenadas dos pontos do arquivo CSV em um quadro de dados e converta-o em SpatialPointsDataFrame
  3. Use a overfunção para determinar em qual polígono ele se encaixa.

Para fazer isso, preciso especificar o proj4stringcarregamento do shapefile na etapa 1 e também transformar as coordenadas do arquivo CSV no mesmo sistema de projeção usando a spTransformfunção antes de aplicar a overfunção na etapa 3, pois exige que os pontos e polígonos sejam necessários. estar sob o mesmo sistema de projeção.

Alguma idéia sobre qual deve ser o valor correto para o conteúdo do arquivo proj mostrado acima?

coordinate-system shapefile r proj Moustafa Alzantot
fonte
Se seu shapefile já possui a projeção definida, use "readOGR" no pacote rgdal. Este pacote é um invólucro para o GDAL e realmente substitui a funcionalidade de leitura / gravação do shapefile no maptools. Esta função lida com todos os tipos de topologia e mantém as informações da projeção.
Jeffrey Evans
Quando tento loadign o arquivo de forma usando readOGRa função I a sempre obter não pode abrir arquivo de erro
Moustafa Alzantot
OK, agora consegui ler o arquivo usando o readOGR. A summaryfunção do SpatialPolygonDataFrameobjeto me forneceu o valor correto para oproj4string
Moustafa Alzantot
Bem, sem detalhes sobre como você está usando a função, é difícil ajudá-lo! Parte da sintaxe é o diretório em que os dados residem e você não precisa da extensão .shp no nome do arquivo. Algo como readOGR (getwd (), "YourShape") deve funcionar se você tiver seu diretório de trabalho definido no mesmo local que seu arquivo shepfile.
Jeffrey Evans
@JeffreyEvans Obrigado, funcionou agora e eu usei-o para obter o proj4string
Moustafa Alzantot

Respostas:

14

O proj4string é uma string PROJ4 crs válida .

consulte Como posso obter a string proj4 ou o código EPSG de um arquivo .prj do shapefile? e Shapefile PRJ para a tabela de pesquisa PostGIS SRID?

em resumo:

  • Você pode usar gdalinfo como na primeira referência ou as ligações GDAL Python como na segunda referência.

Ou

  • vá para Prj2EPSG (um serviço simples para converter informações conhecidas de projeção de texto de arquivos .prj em códigos EPSG padrão)
  • Digite o conteúdo do seu arquivo prj

insira a descrição da imagem aqui

  • o resultado é EPSG: 27700, portanto, uma primeira versão da sequência PROJ4 é

    " + init = epsg: 27700 "

`Ou

insira a descrição da imagem aqui

  • clique no Proj4 e a sequência completa do PROJ4 é:

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

gene
fonte
10

Aqui está um site muito útil para recuperar o código EPSG para uma determinada projeção. No seu caso, a projeção é "EPSG: 27700". Se você tiver projeções definidas para o shapefile, poderá atribuir a projeção ao criar o SpatialPointsDataFrame e, em seguida, usar a definição de projeção do seu shapefile importado. O uso de "readOGR" no pacote rgdal manterá as definições de projeção. Aqui está um exemplo de atribuição e extração de cadeias de coordenadas nos dados da classe 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))
Jeffrey Evans
fonte