Skip to content

Commit

Permalink
Use as(sf::st_crs(wkt), "CRS") #466
Browse files Browse the repository at this point in the history
  • Loading branch information
Jean-Romain committed Sep 3, 2021
1 parent 937d258 commit 2e73434
Show file tree
Hide file tree
Showing 3 changed files with 15 additions and 64 deletions.
67 changes: 11 additions & 56 deletions R/projection.R
Original file line number Diff line number Diff line change
Expand Up @@ -204,7 +204,9 @@ setMethod("projection", "LASheader", function(x, asText = TRUE)
if (use_epsg(x) && epsg(x) != 0L)
proj4 <- epsg2CRS(epsg(x))
else if (use_wktcs(x) && wkt(x) != "")
{
proj4 <- wkt2CRS(wkt(x))
}
else
proj4 <- sp::CRS()

Expand Down Expand Up @@ -425,82 +427,35 @@ use_epsg.LASheader <- function(x) {

epsg2CRS <- function(epsg, fail = FALSE)
{
if (rgdal::new_proj_and_gdal())
{
crs <- tryCatch(
{
sfcrs <- suppressWarnings(sf::st_crs(paste0("EPSG:", epsg)))
return(suppressWarnings(sp::CRS(SRS_string = sfcrs$wkt)))
},
error = function(e)
{
if (!fail)
return(sp::CRS(NA_character_))
else
stop(paste("Invalid epsg code", epsg), call. = FALSE)
})
}
else # nocov start
{
proj <- epsg2proj(epsg, fail)
crs <- sp::CRS(proj)
} # nocov end

return(crs)
}

wkt2CRS <- function(wkt, fail = FALSE)
{
if (!rgdal::new_proj_and_gdal())
{
stop("WKT strings are no longer supported in lidR with old versions of GDAL/PROJ", call. = FALSE)
}

crs <- tryCatch(
{
sp::CRS(SRS_string = wkt)
crs <- suppressWarnings(sf::st_crs(paste0("EPSG:", epsg)))
return(suppressWarnings(as(crs, "CRS")))
},
error = function(e)
{
if (!fail)
return(sp::CRS(NA_character_))
else
stop("Invalid WKT string", call. = FALSE)
stop(paste("Invalid epsg code", epsg), call. = FALSE)
})

return(crs)
}

# nocov start
epsg2proj <- function(epsg, fail = FALSE)
wkt2CRS <- function(wkt, fail = FALSE)
{
tryCatch(
crs <- tryCatch(
{
crs <- sp::CRS(SRS_string = glue::glue("EPSG:{epsg}"))
crs@projargs
as(sf::st_crs(wkt), "CRS")
},
error = function(e)
{
if (!fail)
return(NA_character_)
return(sp::CRS(NA_character_))
else
stop("Invalid epsg code", call. = FALSE)
stop("Invalid WKT string", call. = FALSE)
})
}

wkt2proj <- function(wkt, fail = FALSE)
{
tryCatch(
{
crs <- sp::CRS(SRS_string = glue::glue("EPSG:{epsg}"))
crs@projargs
},
error = function(e)
{
if (!fail)
return(NA_character_)
else
stop("Invalid WKT", call. = FALSE)
})
return(crs)
}
# nocov end
6 changes: 2 additions & 4 deletions tests/testthat/test-las_check.R
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,7 @@ LASfile <- system.file("extdata", "extra_byte.laz", package = "rlas")
las2 <- readLAS(LASfile, select = "xyz")
las2@header@PHB$`Global Encoding`$WKT = TRUE

if(rgdal::new_proj_and_gdal())
wkt(las2) <- "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-81],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"
wkt(las2) <- "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-81],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"

las2@proj4string <- sp::CRS()

Expand Down Expand Up @@ -81,8 +80,7 @@ test_that("las_check CRS specific test", {
las2 <- las0
las2@header@PHB$`Global Encoding`$WKT = TRUE

if(rgdal::new_proj_and_gdal())
wkt(las2) <- "PROJCS[\"RD_New\",GEOGCS[\"GCS_Amersfoort\",DATUM[\"D_Amersfoort\",SPHEROID[\"Bessel_1841\",6377397.155,299.1528128]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Double_Stereographic\"],PARAMETER[\"False_Easting\",155000.0],PARAMETER[\"False_Northing\",463000.0],PARAMETER[\"Central_Meridian\",5.38763888888889],PARAMETER[\"Scale_Factor\",0.9999079],PARAMETER[\"Latitude_Of_Origin\",52.1561605555556],UNIT[\"Meter\",1.0]]"
wkt(las2) <- "PROJCS[\"RD_New\",GEOGCS[\"GCS_Amersfoort\",DATUM[\"D_Amersfoort\",SPHEROID[\"Bessel_1841\",6377397.155,299.1528128]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Double_Stereographic\"],PARAMETER[\"False_Easting\",155000.0],PARAMETER[\"False_Northing\",463000.0],PARAMETER[\"Central_Meridian\",5.38763888888889],PARAMETER[\"Scale_Factor\",0.9999079],PARAMETER[\"Latitude_Of_Origin\",52.1561605555556],UNIT[\"Meter\",1.0]]"

las2@header@PHB$`Global Encoding`$WKT <- FALSE

Expand Down
6 changes: 2 additions & 4 deletions tests/testthat/test-projection.R
Original file line number Diff line number Diff line change
Expand Up @@ -89,10 +89,7 @@ test_that("epsg<- works", {
test_that("Set an invalid code or WKT fails", {
expect_error({ epsg(las) <- 200800 }, "Invalid epsg")

if (rgdal::new_proj_and_gdal())
expect_error({ wkt(las) <- "INVALID" }, "Invalid WKT")
else
expect_error({ wkt(las) <- "INVALID" }, "WKT strings are no longer supported in lidR with old versions of GDAL/PROJ")
expect_error({ wkt(las) <- "INVALID" }, "Invalid WKT")
})

test_that("#323 do not segfault", {
Expand All @@ -101,3 +98,4 @@ test_that("#323 do not segfault", {
wkt = "PROJCS[\"NAD83 (2011) / Conus Albers\",GEOGCS[\"GCS_NAD_1983_2011\",DATUM[\"D_NAD_1983_2011\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101,AUTHORITY[\"EPSG\",7019]],AUTHORITY[\"EPSG\",1116]],PRIMEM[\"Greenwich\",0.0,AUTHORITY[\"EPSG\",8901]],UNIT[\"Degree\",0.0174532925199433,AUTHORITY[\"EPSG\",9102]],AUTHORITY[\"EPSG\",6318]],PROJECTION[\"Albers\",AUTHORITY[\"Esri\",43007]],PARAMETER[\"False_Easting\",0.0,AUTHORITY[\"Esri\",100001]],PARAMETER[\"False_Northing\",0.0,AUTHORITY[\"Esri\",100002]],PARAMETER[\"Central_Meridian\",-96.0,AUTHORITY[\"Esri\",100010]],PARAMETER[\"Standard_Parallel_1\",29.5,AUTHORITY[\"Esri\",100025]],PARAMETER[\"Standard_Parallel_2\",45.5,AUTHORITY[\"Esri\",100026]],PARAMETER[\"Latitude_Of_Origin\",23.0,AUTHORITY[\"Esri\",100021]],UNIT[\"Meter\",1.0,AUTHORITY[\"EPSG\",9001]]],VERTCS[\"NAVD_1988\",VDATUM[\"North_American_Vertical_Datum_1988\",AUTHORITY[\"EPSG\",5103]],PARAMETER[\"Vertical_Shift\",0.0,AUTHORITY[\"Esri\",100006]],PARAMETER[\"Direction\",1.0,AUTHORITY[\"Esri\",100007]],UNIT[\"Meter\",1.0,AUTHORITY[\"EPSG\",9001]],AUTHORITY[\"EPSG\",5703]]"
expect_error(lidR:::wkt2CRS(wkt), NA)
})

0 comments on commit 2e73434

Please sign in to comment.