Koda
Koda

Reputation: 177

R Projection issue with raster and polygon not overlapping even though they do

I need to crop a raster (Sentinel 2 image) by a polygon (USGS HUC data). The raster's crs is: Deprecated Proj.4 representation: +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs . The polygon's crs is "GEOGCRS["NAD83",\n DATUM["North American Datum 1983",\n ..."

I have tried spatial transforming the polygon to the raster's crs and also tried projecting the raster in the polygon's crs. The extent of both items overlap but are in differing units (UTM vs lat/long).

What is the trick to changing the projections so that the two overlap?

Here is the code:

UpCarson <- read_sf("UpperCarson.shp")
crs(UpCarson)
> crs(Mar17_project)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=NAD83 +no_defs 
WKT2 2019 representation:
GEOGCRS["unknown",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6269]],
    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]]]] 

Mar17 <- raster("S2A2A_20230317_070_sen2r_NDSI_10.tif")
crs(Mar17)
> crs(Mar17)
Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
WKT2 2019 representation:
PROJCRS["unknown",
    BASEGEOGCRS["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]]],
    CONVERSION["UTM zone 10N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]],
        ID["EPSG",16010]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]] 

Mar17_project <- projectRaster(Mar17, crs = crs(UpCarson))
crs(Mar17_project)

> crs(Mar17_project)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +datum=NAD83 +no_defs 
WKT2 2019 representation:
GEOGCRS["unknown",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6269]],
    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]]]] 

##Crop the raster by the polygon
Mar17_crop <- crop(Mar17_project, UpCarson)
> Mar17_crop <- crop(Mar17_project, UpCarson)
Error in .local(x, y, ...) : extents do not overlap

Upvotes: 0

Views: 158

Answers (1)

Koda
Koda

Reputation: 177

The polygon was brought in as a sf but the spatial transformation requires that it is a sp. Converting the sf to a sp allowed for the spatial transformation to take place.

UpCarson_sp <- as(UpCarson, "Spatial")
##crs taken directly from the raster crs
UpCarson_project <- spTransform(UpCarson, CRSobj = "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs")

Upvotes: 0

Related Questions