sp :: proj4string эквивалент в структуре PROJ6?

Question

Меня смущает переход от PROJ4 к PROJ6 в (относительно) последних версиях пространственных пакетов R (и gdal в этом отношении). Раньше я использовал, sp::proj4string()чтобы получить код проекции / аргумент CRS от spобъекта и использовать этот код в своих сценариях. Теперь я получаю предупреждение (в котором говорится, что мой прогноз может быть неточным ). Я явно что-то не так делаю. Какую функцию следует использовать для получения кода проекции без предупреждения? Я лично предпочел бы код EPSG, если бы такой уровень простоты был возможен ...

library(sp)

x <- SpatialPolygons(
            list(Polygons(
                 list(Polygon(cbind(c(1,2,3,4,1), c(3,8,10,12,3)))), 
                 ID = 1)), 
            proj4string = CRS(SRS_string = "EPSG:4326"))

proj4string(x)
#> Warning in proj4string(x): CRS object has comment, which is lost in output
#> [1] "+proj=longlat +datum=WGS84 +no_defs"

^{Создано 18 сентября 2020 г. пакетом REPEX (v0.3.0)}

См. Также ссылки и комментарии в этом вопросе.

Answer 1

Я также был весьма озадачен отказом от PROJ Strings и тем, что использовать вместо этого, пока не посмотрел это видео с выступления, проведенного в прошлом году на FOSS4G.

Вот мои основные выводы:

Почему PROJ и GDAL отходят от PROJ-строк

В PROJ4 каждое преобразование выполняется с использованием «концентрационного подхода». Так что сначала все будет преобразовано в WGS84, а только потом преобразовано в целевую проекцию.

Такой подход может привести к большим ошибкам в процессе преобразования, чем прямое преобразование между двумя проекциями. В худшем случае необходимое преобразование даже невозможно при таком подходе. Поскольку строки PROJ могут быть неоднозначными, прямые преобразования часто невозможны.

Пример из разговора:

И старые австралийские данные GDA94, и новые GDA2020 основаны на WGS84. Однако один из них использует опорный кадр ITRF2014 в Эпоху 2020.0, а другой ITRF92 в Эпоху 1994.0. Эти различия не могут быть зафиксированы в строке PROJ (насколько я понял), и поэтому все проекции имеют почти ту же строку PROJ, что и WGS84. Преобразование между обоими с использованием подхода ступицы не приведет к изменению координат. Хотя на самом деле Австралия переместилась на 1,8 метра между 1994 и 2020 годами, и преобразование датума должно это отразить.

При прямом преобразовании таких ошибок не будет.

Что вам следует использовать вместо

Вместо этого вы должны в идеале использовать однозначные определения для прогнозов, такие как коды EPSG или определения в формате WKT2.

Однако, если вы не собираетесь выполнять преобразование датума (или не заботитесь о погрешности в несколько метров или более), вы все равно можете использовать PROJ-строки, но имейте в виду, что лучше всего использовать коды EPSG или WKT2.

Как это реализовано в sp

CRS-объекты в sp теперь имеют дополнительное поле, commentв котором сохраняется WKT2-представление CRS.

Итак, рабочий процесс с лучшими практиками spсейчас выглядит следующим образом (взято и немного расширено отсюда ):

# Define the CRS using an EPSG Code
x <- CRS(SRS_string='EPSG:4326')

# Display the stored CRS using comment()
cat(comment(x), "\n")

# Store the wkt in a variable
wkt <- comment(x)

# Use this to assign the CRS of another sp-object
y <- CRS(SRS_string = wkt)

В общем, sfимеет гораздо более удобные функции для доступа и обработки crs, чем sp. Вы можете узнать больше о том, как sfэто работает, по той же ссылке, указанной выше.

Answer 2

Proj4strings постепенно сокращается. Благодаря пакету sf и его разработчикам есть простой трюк для преобразования proj4string в формат WKT.

Просто игнорируйте предупреждения о данных, исходящие от rgdal, и, по сути, вы можете просто отключить их, потому что они действительно надоедают. Либо добавьте options("rgdal_show_exportToProj4_warnings"="none")в файл Rprofile.site, либо выдайте его перед добавлением каких-либо библиотек.

Например, с географической проекцией WGS84 использование sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")приведет к

Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +no_defs 
  wkt:
GEOGCRS["unknown",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]

Из блога R-Space R Space следует за разработкой GDAL и PROJ, и разработчики утверждают:

Использование так называемых PROJ4-строк (например, + proj = longlat + datum = WGS84) не рекомендуется, они больше не предлагают достаточного описания систем координат; использование + init = epsg: XXXX приводит к предупреждению

Answer 3

Это предупреждение никак не влияет на результат. Как предлагали другие, пакет sf может быть полезен для более чистого и простого кодирования EPSG. Вот ваш адаптированный код:

library(sf)
library(mapview)
coords <- list(rbind(c(1,3), c(2,8), c(3,10), c(4,12), c(1,3)))
x2 <- st_sf(st_geometry(st_polygon(x = coords)), crs="EPSG:4326")
x2$ID <- 1
mapview(x2)