DOI-USGS / meddle

Tools for metadata creation and data releases
Other
0 stars 10 forks source link

mismatched CRS in get_states() #84

Closed aappling-usgs closed 5 months ago

aappling-usgs commented 1 year ago

Same function as issue #76

meddle:::get_states()
Error: arguments have different crs

> traceback()
9. stop("arguments have different crs", call. = FALSE)
8. FUN(X[[i]], ...)
7. vapply(dots[-1L], function(x) {
     if (st_crs(x) != crs0)
     stop("arguments have different crs", call. = FALSE)
     TRUE ...
6. chk_equal_crs(dots)
5. rbind(deparse.level, ...)
4. rbind(us_48, us_ak) at extract_spatial.R#255
3. rbind(rbind(us_48, us_ak), us_hi) at extract_spatial.R#255
2. rbind(rbind(rbind(us_48, us_ak), us_hi), us_pr) at extract_spatial.R#255
1. meddle:::get_states()

With debug(meddle:::get_states); meddle:::get_states(), I find that before the setting of st_crs for us_pr, we have these Coordinate Reference Systems set (calling sf::st_crs(us_area) on each area:

Browse[2]> sf::st_crs(us_pr)
Coordinate Reference System: NA

Browse[2]> sf::st_crs(us_48)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Browse[2]> sf::st_crs(us_ak)
Coordinate Reference System:
  User input: +proj=longlat +ellps=clrk66 +no_defs +type=crs 
  wkt:
GEOGCRS["unknown",
    DATUM["Unknown based on Clarke 1866 ellipsoid",
        ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    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]]]]

Browse[2]> sf::st_crs(us_hi)
Coordinate Reference System:
  User input: +proj=longlat +ellps=clrk66 +no_defs +type=crs 
  wkt:
GEOGCRS["unknown",
    DATUM["Unknown based on Clarke 1866 ellipsoid",
        ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    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]]]]
aappling-usgs commented 1 year ago

I added these two lines within the get_states() function (before rbind()) and it seems to work, but I haven't yet figured out how to submit a PR on github.com/DOI-USGS.

us_hi <- sf::st_transform(us_hi, crs = 4326)
us_ak <- sf::st_transform(us_ak, crs = 4326)