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

ViSa

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()

enter image description here

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()

enter image description here

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)

enter image description here

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]]
dww

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))

Collected from the Internet

Please contact [email protected] to delete if infringement.

edited at
0

Comments

0 comments
Login to comment

Related

How to make future projections without data ? (R)

shapefile to a raster .tif file R

R: Aggregating raster to shapefile polygons

R and Export raster extent as esri shapefile

convert shapefile of multiple polygons to raster in R

Overlap (or convert) a raster in shapefile with R software

Overlaying shapefile over multiple raster plots in R

How to make data frame into raster object R

calculating area of most suitable raster habitat in R

How to make projections with future dates using Redshift

How to make LINQ-to-Objects handle projections?

How to make different orthogonal projections Pyglet?

How do I get this raster and this shapefile on the same projection?

How to write a loop for creating cropped raster for every id of a shapefile with a raster base?

How do I use rasterio/python to mask a raster using a shapefile, to set the raster pixels inside the polygons to zero?

How to extract lists of raster into separate raster in R

masking big raster using shapefile

How to make smooth circles on basemap projections in Matplotlib by Python

How to rotate an image R raster

How to extrapolate a raster using in R

How to solve problem with reading a shapefile in R

R Shiny, how to display a shapefile (map) in a tab?

How to limit boundaries when plotting a shapefile in R

How to use CAST package for a shapefile (polygons) in R?

Finding differences in projections in R and projections in ArcMap

How to make my HTML suitable for Mobile View

How to make raster from unique values?

How to make this "cropping"/background color changing function more efficient?

How to make pandas show the entire dataframe without cropping it by columns?