Rodrigo-NH / GeoJsonIsRight

Test to show how right/wrong is GeoJson about having Logintude first
BSD 3-Clause "New" or "Revised" License
1 stars 0 forks source link

R version (for mutual learning) #1

Open mdsumner opened 2 years ago

mdsumner commented 2 years ago

thanks for this cool example!

Here's a stab at an R version, not sure my logic is sound yet but this is great for me to compare R and python approaches

library(dbplyr)
library(dplyr)

##dbfile <- libproj::libproj_configuration()$db_path
dbfile <- "/usr/local/lib/R/site-library/libproj/proj/proj.db"
db <- DBI::dbConnect(RSQLite::SQLite(),  dbfile)
epsg <- tbl(db, sql("SELECT code FROM crs_view WHERE auth_name = 'EPSG'")) |> 
  collect() |> 
  mutate(auth = sprintf("EPSG:%i", code))

compare_axis <- \(.x) {
  wkt_lines <- vapour::vapour_srs_wkt(.x)
  #if (grepl("AXIS.*OTHER", wkt_lines)) return(NA)
  lonax <- c(attr(gregexpr("AXIS.*EAST", wkt_lines)[[1]], "match.length"), 
             attr(gregexpr("AXIS.*LONGITUDE", wkt_lines)[[1]], "match.length"))
  latax <- c(attr(gregexpr("AXIS.*NORTH", wkt_lines)[[1]], "match.length"), 
             attr(gregexpr("AXIS.*LATITUDE", wkt_lines)[[1]], "match.length"))
  if (all(latax < 1) || all(lonax < 1)) return(NA)
  min(lonax[lonax > 0]) < min(latax[latax > 0])
}

epsg <- epsg |> mutate(lonlat =  purrr::map_lgl(epsg$auth, compare_axis))
epsg
#> # A tibble: 7,049 × 3
#>     code auth      lonlat
#>    <int> <chr>     <lgl> 
#>  1  3819 EPSG:3819 FALSE 
#>  2  3821 EPSG:3821 FALSE 
#>  3  3822 EPSG:3822 NA    
#>  4  3823 EPSG:3823 NA    
#>  5  3824 EPSG:3824 FALSE 
#>  6  3887 EPSG:3887 NA    
#>  7  3888 EPSG:3888 NA    
#>  8  3889 EPSG:3889 FALSE 
#>  9  3906 EPSG:3906 FALSE 
#> 10  4000 EPSG:4000 NA    
#> # … with 7,039 more rows

epsg |> filter(!is.na(lonlat)) |> summarize(mean(lonlat))
#> # A tibble: 1 × 1
#>   `mean(lonlat)`
#>            <dbl>
#> 1          0.638

system("proj")
# Rel. 8.2.0, November 1st, 2021

Created on 2022-08-16 by the reprex package (v2.0.1)

Rodrigo-NH commented 2 years ago

Thanks! I'm still 'R' illiterate and your example is (yet another) incentive to get myself into the language. What most annoyed me is the diversity in the nomenclature for the same thing (NORTH/NORTHING/LATITUDE and EAST/EASTING/LONGITUDE), perhaps subtle details I'm not aware, considering EPSG was some coordinated initiative.

Other thing that's great about doing it is knowing more what's inside proj.db. Interesting tables.

[('alias_name',), ('authority_to_authority_preference',), ('axis',), ('celestial_body',), ('compound_crs',), ('concatenated_operation',), ('concatenated_operation_step',), ('conversion_method',), ('conversion_param',), ('conversion_table',), ('coordinate_operation_method',), ('coordinate_system',), ('deprecation',), ('ellipsoid',), ('extent',), ('geodetic_crs',), ('geodetic_datum',), ('geodetic_datum_ensemble_member',), ('geoid_model',), ('grid_alternatives',), ('grid_packages',), ('grid_transformation',), ('helmert_transformation_table',), ('metadata',), ('other_transformation',), ('prime_meridian',), ('projected_crs',), ('scope',), ('sqlite_stat1',), ('sqlite_stat4',), ('supersession',), ('unit_of_measure',), ('usage',), ('versioned_auth_name_mapping',), ('vertical_crs',), ('vertical_datum',), ('vertical_datum_ensemble_member',)]

mdsumner commented 2 years ago

east and north, easting and northing, is a misnomer - x and y makes sense but east and north do not generally