How to adjust raster & shapefile projections in r to make it suitable for cropping?

Question

I am new to geospatial data & trying to crop a tif file raster object using a shapefile by referring https://www.youtube.com/watch?v=UP2Za1TizOc.

I have tried below code by referring the above video & seems like there is a projection issue:

library(tidyverse)
library(sf)
library(stars)

ind_outline <- sf::st_read("path\\polymap15m_area.shp")
ind_outline


Simple feature collection with 314 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2815341 ymin: 2177499 xmax: 5678865 ymax: 5444567
Projected CRS: LCC_WGS84
First 10 features:
   Id Line_Width                       geometry
1   0       1875 POLYGON ((5547296 2230982, ...
2   0       1875 POLYGON ((5560180 2232030, ...
3   0       1875 POLYGON ((5549993 2253154, ...
4   0       1875 POLYGON ((5542651 2256150, ...
5   0       1875 POLYGON ((5533962 2260494, ...
6   0       1875 POLYGON ((5523175 2264240, ...
7   0       1875 POLYGON ((3223295 2294948, ...
8   0       1875 POLYGON ((5502051 2315325, ...
9   0       1875 POLYGON ((5522126 2328209, ...
10  0       1875 POLYGON ((5480027 2338995, ...

shapefile link: https://github.com/johnsnow09/covid19-df_stack-code/blob/main/polymap15m_area.shp

ind_outline %>% 
        st_as_sf() %>% 
        ggplot() +
        geom_sf()

ind_region_stars <- stars::read_stars("../gt30e060n40.tif")
ind_region_stars

stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                 Min. 1st Qu. Median     Mean 3rd Qu. Max.
gt30e060n40.tif   130     793   1052 1302.186    1648 4795
dimension(s):
  from   to offset       delta refsys point values x/y
x    1 4800     60  0.00833333 WGS 84 FALSE   NULL [x]
y    1 6000     40 -0.00833333 WGS 84 FALSE   NULL [y]

ggplot() +
        geom_stars(data = ind_region_stars) +
        scale_fill_viridis_c()

Issue: When I try to overlay one on other then it doesn't work

ggplot() +
        geom_stars(data = ind_region_stars) +
        scale_fill_viridis_c() +
        geom_sf(data = ind_outline, alpha = 0)

Error in cropping:

ind_region_stars_cropped <- st_crop(ind_region_stars, ind_outline)

Error in st_crop.stars(ind_region_stars, ind_outline) : for cropping, the CRS of both objects have to be identical

UPDATE:

st_crs(ind_outline)

Coordinate Reference System:
  User input: LCC_WGS84 
  wkt:
PROJCRS["LCC_WGS84",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",24,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",80,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",12.472944,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",35.172806,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",4000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",4000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    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]]]]

st_crs(ind_region_stars)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    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]],
    ID["EPSG",4326]]

you simply need to reproject both into the same desired CRS. — dww, Jun 29 '21 at 10:50
yes but I was not able to figure out lcc .. I guess this will work: https://stackoverflow.com/questions/30287065/convert-lambert-conformal-conic-projection-to-wgs84-in-r — ViSa, Jun 29 '21 at 10:51
if you want to set polygons to same crs as raster then do that explicitly - ie reproject to `ind_outline <- st_transform(ind_outline, crs = st_crs(ind_region_stars))` — dww, Jun 29 '21 at 10:55
@dww I have tried below code but getting error; `crs <- CRS("+proj=lcc +lat_1=30 +lat_2=60 +lat_0=38 +lon_0=126 +datum=WGS84")` `ind_outline_crs <- SpatialPoints(ind_outline, proj4string=crs)` — ViSa, Jun 29 '21 at 10:55
"getting error" not very useful. Please read https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example — dww, Jun 29 '21 at 10:57
thanks alot @dww your code worked. If you would add this in the answer then I would accept it & close this post. — ViSa, Jun 29 '21 at 10:58

score 1 · Accepted Answer · answered Jun 29 '21 at 11:02

1

To set polygons to same crs as raster then you should do that explicitly using

ind_outline <- st_transform(ind_outline, crs = st_crs(ind_region_stars))

answered Jun 29 '21 at 11:02

dww

30,425
5
68
111

How to adjust raster & shapefile projections in r to make it suitable for cropping?

1 Answers1