tests/testthat/test-gridmapping-prj.R

In hypertidy/ncmeta: Straightforward 'NetCDF' Metadata

context("prj and grid mappings")

test_that("nc_grid_mapping_atts", {
  
  nc <- system.file("extdata/S2008001.L3m_DAY_CHL_chlor_a_9km.nc", package = "ncmeta")
  expect_warning(gm <- nc_grid_mapping_atts(nc),
                 paste("No variables with a grid mapping found.\n",
                       "Defaulting to WGS84 Lon/Lat"))
  expect_equal(stats::setNames(gm$value, gm$name), list(grid_mapping_name = "latitude_longitude",
                        semi_major_axis = 6378137,
                        inverse_flattening = 298.257223563,
                        longitude_of_prime_meridian = 0))
  
  nc <- system.file("extdata/daymet_sample.nc", package = "ncmeta")
  gm <- nc_grid_mapping_atts(nc)
  
  expect_true(all(list(grid_mapping_name = "lambert_conformal_conic",
                  longitude_of_central_meridian = -100,
                  latitude_of_projection_origin = 42.5,
                  false_easting = 0,
                  false_northing = 0,
                  standard_parallel = c(25, 60),
                  semi_major_axis = 6378137,
                  inverse_flattening = 298.257223563,
                  longitude_of_prime_meridian = 0) %in% stats::setNames(gm$value, gm$name)))
  
  expect_is(nc_grid_mapping_atts(ncmeta::nc_atts(nc)), "data.frame")
  
  expect_is(nc_grid_mapping_atts(RNetCDF::open.nc(nc)), "data.frame")
  
  gm2 <- nc_grid_mapping_atts(nc, data_variable = "prcp")
  
  expect_equal(nrow(gm), nrow(gm2))
  
  expect_warning(nc_grid_mapping_atts(nc, data_variable = "borked"),
                 "no grid_mapping attribute found for this variable")
})

test_that("nc_prj_to_gridmapping returns an empty list if no mapping exists", {
  p <- ""

  expect_warning(crs <- nc_prj_to_gridmapping(p), 
                 "not a valid crs, returning an empty tibble")

  expect_equal(names(crs), c("name", "value"))
  expect_equal(nrow(crs), 0)
})

test_that("wgs 84 lat lon", {
  p <- "+proj=longlat +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name="latitude_longitude",
              longitude_of_prime_meridian = 0,
              semi_major_axis = 6378137,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)

  crs <- stats::setNames(crs$value, crs$name)
  
  expect_equal(crs, c[names(crs)])
})

test_that("NAD27 lat lon", {
  p <- "+proj=longlat +datum=NAD27 +no_defs"
  
  p2 <- "+proj=longlat +a=6378206.4 +f=0.00339007530392762 +pm=0 +no_defs"
  
  c <- list(grid_mapping_name="latitude_longitude",
            longitude_of_prime_meridian = 0,
            semi_major_axis = 6378206.4,
            inverse_flattening = 294.978698214)
  
  prj <- nc_gm_to_prj(c)
  
  expect_equal(prj, p2)
  
  crs <- nc_prj_to_gridmapping(p)
  
  crs <- stats::setNames(crs$value, crs$name)
  
  expect_equal(crs, c[names(crs)])
})

test_that("albers equal area epsg:5070", {
  p <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +units=m +lat_0=23 +lon_0=-96 +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "albers_conical_equal_area",
              longitude_of_central_meridian = -96,
              latitude_of_projection_origin = 23,
              false_easting = 0.0,
              false_northing = 0.0,
              standard_parallel = c(29.5, 45.5),
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563,
              longitude_of_prime_meridian = 0)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)
  
  expect_equal(crs, c[names(crs)])
})

test_that("albers equal area epsg:5070 with datum instead of a b", {
  p <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +units=m +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

  c <- list(grid_mapping_name = "albers_conical_equal_area",
            longitude_of_central_meridian = -96,
            latitude_of_projection_origin = 23,
            false_easting = 0.0,
            false_northing = 0.0,
            standard_parallel = c(29.5, 45.5),
            semi_major_axis = 6378137.0,
            inverse_flattening = 298.257222101,
            longitude_of_prime_meridian = 0)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs, c[names(crs)])
})

test_that("Azimuthal Equidistant", {
  p <- "+proj=aeqd +lat_0=30 +lon_0=-40 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "azimuthal_equidistant",
              longitude_of_projection_origin = -40,
              latitude_of_projection_origin = 30,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)

})

test_that("lambert conformal conic daymet", {
  p <- "+proj=lcc +lat_1=25 +lat_2=60 +x_0=0 +y_0=0 +units=m +lat_0=42.5 +lon_0=-100 +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "lambert_conformal_conic",
              longitude_of_central_meridian = -100.0,
              latitude_of_projection_origin = 42.5,
              false_easting = 0.0,
              false_northing = 0.0,
              standard_parallel = c(25.0, 60.0),
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)

})

test_that("lambert_azimuthal_equal_area", {
  p <- "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +units=m +a=6371228 +b=6371228 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "lambert_azimuthal_equal_area",
              longitude_of_projection_origin = 0,
              latitude_of_projection_origin = 90,
              false_easting = 0.0,
              false_northing = 0.0,
              semi_major_axis = 6371228,
              semi_minor_axis = 6371228,
              longitude_of_prime_meridian = 0.0)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  c <- list(grid_mapping_name = "lambert_azimuthal_equal_area",
              longitude_of_projection_origin = 0,
              latitude_of_projection_origin = 90,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6371228)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)

})


test_that("lambert_cylindrical_equal_area", {
  p <- "+proj=cea +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "lambert_cylindrical_equal_area",
              longitude_of_central_meridian = 0,
              standard_parallel=0,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)

})

test_that("mercator", {
  p <- "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "mercator",
              longitude_of_projection_origin = 0,
              standard_parallel=0,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)

  p <- "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "mercator",
              longitude_of_projection_origin = 0,
              scale_factor_at_projection_origin=1,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)


})

test_that("oblique_mercator", {
  p <- "+proj=omerc +lat_0=45.3091666666667 +lonc=-86 +k=0.9996 +alpha=337.25556 +gamma=337.25556 +no_uoff +x_0=2546731.496 +y_0=-4354009.816 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "oblique_mercator",
              azimuth_of_central_line = 337.25556,
              longitude_of_projection_origin = -86,
              latitude_of_projection_origin = 45.3091666666667,
              scale_factor_at_projection_origin = 0.9996,
              false_easting = 2546731.496,
              false_northing = -4354009.816,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)

})

test_that("orthographic", {
  p <- "+proj=ortho +lat_0=30 +lon_0=-40 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "orthographic",
              longitude_of_projection_origin = -40,
              latitude_of_projection_origin = 30,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)

})

test_that("polar_stereographic", {
  p <- "+proj=stere +lat_0=-90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "polar_stereographic",
              straight_vertical_longitude_from_pole = 0,
              latitude_of_projection_origin = -90,
              scale_factor_at_projection_origin = 1,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  # crs <- nc_prj_to_gridmapping(p)
  # crs <- stats::setNames(crs$value, crs$name)
  #
  # expect_equal(crs[names(c)], c)

  p <- "+proj=stere +lat_0=-90 +lon_0=0 +lat_ts=-71 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "polar_stereographic",
              straight_vertical_longitude_from_pole = 0,
              latitude_of_projection_origin = -90,
              standard_parallel = -71,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  # crs <- nc_prj_to_gridmapping(p)
  # crs <- stats::setNames(crs$value, crs$name)
  #
  # expect_equal(crs[names(c)], c)

})

test_that("sinusoidal", {
  p <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "sinusoidal",
              longitude_of_projection_origin = 0,
              false_easting = 0.0,
              false_northing = 0.0,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

})

test_that("stereographic", {
  p <- "+proj=stere +lat_0=46.5 +lon_0=-66.5 +k=0.999912 +x_0=2500000 +y_0=7500000 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"

  c <- list(grid_mapping_name = "stereographic",
              longitude_of_projection_origin = -66.5,
              latitude_of_projection_origin = 46.5,
              scale_factor_at_projection_origin = 0.999912,
              false_easting = 2500000,
              false_northing = 7500000,
              longitude_of_prime_meridian = 0.0,
              semi_major_axis = 6378137.0,
              inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)


})

test_that("transverse_mercator", {
  p <- "+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=0.99982 +x_0=600000 +y_0=750000 +units=m +a=6378137 +f=0.00335281066474748 +pm=0 +no_defs"
  
  c <- list(grid_mapping_name = "transverse_mercator",
            longitude_of_central_meridian = -8,
            latitude_of_projection_origin = 53.5,
            scale_factor_at_central_meridian = 0.99982,
            false_easting = 600000,
            false_northing = 750000,
            longitude_of_prime_meridian = 0.0,
            semi_major_axis = 6378137.0,
            inverse_flattening = 298.257223563)

  prj <- nc_gm_to_prj(c)

  expect_equal(prj, p)

  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)

  expect_equal(crs[names(c)], c)


})

test_that("spherical", {
  
  p <- "+proj=lcc +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +units=m +lat_0=40.0000076294 +lon_0=-97 +a=6370000 +b=6370000 +pm=0 +no_defs"
  
  # From the National Water Model forcings
  cl <- list(transform_name = "lambert_conformal_conic", 
            grid_mapping_name = "lambert_conformal_conic", 
            standard_parallel = c(30, 60), 
            longitude_of_central_meridian = -97, 
            latitude_of_projection_origin = 40.0000076294, 
            false_easting = 0, 
            false_northing = 0, 
            earth_radius = 6370000)

  prj <- suppressWarnings(nc_gm_to_prj(cl))
  
  expect_equal(prj, p)
  
  crs <- nc_prj_to_gridmapping(p)
  crs <- stats::setNames(crs$value, crs$name)
  
  expect_equal(crs, 
               list(grid_mapping_name = "lambert_conformal_conic", 
                    standard_parallel = c(30, 60), 
                    false_easting = 0, 
                    false_northing = 0, 
                    latitude_of_projection_origin = 40.0000076294, 
                    longitude_of_central_meridian = -97, 
                    semi_major_axis = 6370000, 
                    semi_minor_axis = 6370000,
                    longitude_of_prime_meridian = 0))
 
})

test_that("more spherical", {
  p <- "+proj=lcc +lat_1=30 +lat_2=65 +x_0=-6000 +y_0=-6000 +units=m +lat_0=48 +lon_0=9.75 +a=6371229 +b=6371229 +pm=0 +no_defs"
  
  # from "https://github.com/mdsumner//rasterwise/extdata/EURO-CORDEX_81_DOMAIN000_54/EURO-CORDEX_81_DOMAIN000.nc"
  cl <- list(proj4_params = "+proj=lcc +lat_1=30.00 +lat_2=65.00 +lat_0=48.00 +lon_0=9.75 +x_0=-6000. +y_0=-6000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs", 
             grid_mapping_name = "lambert_conformal_conic", 
             standard_parallel = c(30, 65), 
             longitude_of_central_meridian = 9.75, 
             latitude_of_projection_origin = 48, 
             semi_major_axis = 6371229, 
             inverse_flattening = 0, # This is the BS that caused an Inf!
             false_easting = -6000, 
             false_northing = -6000, 
             `_CoordinateTransformType` = "Projection", 
             `_CoordinateAxisTypes` = "GeoX GeoY")
  

   prj <- suppressWarnings(nc_gm_to_prj(cl))
   
   expect_equal(prj, p)
   
   crs <- nc_prj_to_gridmapping(p)
   crs <- stats::setNames(crs$value, crs$name)
  
   expect_equal(crs, 
                list(grid_mapping_name = "lambert_conformal_conic", 
                     standard_parallel = c(30, 65), 
                     false_easting = -6000, 
                     false_northing = -6000, 
                     latitude_of_projection_origin = 48, 
                     longitude_of_central_meridian = 9.75, 
                     semi_major_axis = 6371229, 
                     semi_minor_axis = 6371229, 
                     longitude_of_prime_meridian = 0))
})