Как получить код EPSG из строки PROJ4, используя R? - PullRequest
Новая
       97

Как получить код EPSG из строки PROJ4, используя R?

1 голос
/ 08 января 2020

Как получить код EPSG из строки PROJ4, используя R ?

У меня есть набор растровых данных (случайным образом распределенных по всему миру, в северном и южном полушарии), с система координат UTM WGS84. Теперь я хотел бы узнать конкретный c код EPSG для каждого набора данных, использующего R. Есть ли у вас какие-либо идеи? 

# Install and load package raster
install.packages('raster')
library(raster)

# Load a specific raster dataset and query the crs
raster_dat <- raster('D:/UnknownUser/GlobalRasterData/raster_dat1')
crs(raster_dat)
# CRS arguments:
#  +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

Может быть, возможно создать код EPSG в динамическом c Кстати, поскольку кажется, что код EPSG всегда является комбинацией <326> для северного полушария и указанного c номера зоны, например, <32> (для Германии), что приводит к коду EPSG <326> + <32> = <32632. 327> + <32> = <32732>.

Теперь я не уверен, верна ли эта теория ( требуется проверка, включая источник, пожалуйста, ) и есть ли более простой способ получить код EPSG из строки PROJ4.

Полезные ссылки, которые я нашел:

URL: https://spatialreference.org/ref/epsg/?page=9&search=wgs+84+utm

URL: https://epsg.io/32632

URL: https://de.wikipedia.org/wiki/UTM-Koordinatensystem# / media / Datei: Utm-zone.jpg

РЕДАКТИРОВАТЬ: 

# Install and load package rgdal
install.packages('rgdal')
library(rgdal)

# following function could be helpful; it creates a dataframe with 3columns; epsg-code, annotation and proj4-string - unfortunately it seems like there is not an epsg code for each utm wgs84 zone...
epsg_list <- rgdal::make_EPSG()

Заранее большое спасибо, ExploreR

1 Ответ

1 голос
/ 09 января 2020

Пакет rgdal также имеет функцию showEPSG(). 


## example data from sf package
nc <- sf::read_sf(system.file("shape/nc.shp", package="sf"))
nc_crs <- raster::crs(nc)
nc_crs
#> [1] "+proj=longlat +datum=NAD27 +no_defs"

rgdal::showEPSG(nc_crs)
#> [1] "4267"

## Double check it with sf::st_crs
sf::st_crs(nc)
#> Coordinate Reference System:
#>   EPSG: 4267 
#>   proj4string: "+proj=longlat +datum=NAD27 +no_defs"

## Your proj4 string
p4_example <- '+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0;'
rgdal::showEPSG(p4_example)
#> [1] "32613"

Создан в 2020-01-09 пакетом Представить (v0.3.0)

...