Error in geos_op2_geom("intersection", x, y, ...) : st_crs(x) == st_crs(y) is not TRUE

[moniqa]

Hi.
I am trying to run my script:

chm_p2r = raster("chm_p2r/rasterize_canopy.vrt")
chm_p2r = readAll(chm_p2r)
ttops_chm_p2r = st_read("ttops.shp")
ttops_chm_p2r = st_cast(ttops_chm_p2r, to = "POINT")
st_crs(ttops_chm_p2r) = 2180
crowns = dalponte2016(chm_p2r, ttops_chm_p2r)()

And then I have an error:

Error in geos_op2_geom("intersection", x, y, ...) : 
  st_crs(x) == st_crs(y) is not TRUE

Maybe chm_p2r does not have 2180 crs and I know that chm_p2r is raster layer and when I try:

hm_p2r_new = st_set_crs(chm_p2r, 2180)

I have:

Error in UseMethod("st_crs<-") : no applicable method for 'st_crs<-' applied to an object of class "c('RasterLayer', 'Raster', 'BasicRaster')

How to set crs to RasterLayer? Or how to solve this problem, because maybe there is some other option?

Thanks in advance!

[technocrat]

Without a reprex (see the FAQ), I haven't tested this, but it should work.

To conform the Coordinate Reference System (CRS) of a raster layer to another raster layer in R, you can use the projectRaster() function from the raster package[1][2]. Here's a step-by-step guide on how to do this:

1. Load the raster package:

library(raster)

2. Load or create the two raster objects, for example, raster1 and raster2. You want to conform the CRS of raster2 to match the CRS of raster1.

3. Use the crs() function to obtain the CRS of raster1:

raster1_crs <- crs(raster1)

4. Use the projectRaster() function to reproject raster2 to match the CRS of raster1:

raster2_conformed <- projectRaster(raster2, crs = raster1_crs)

Now, raster2_conformed has the same CRS as raster1, and you can perform further analysis or visualization with these raster layers[1][2].

Citations:
[1] Introduction to Geospatial Raster and Vector Data with R: Reproject Raster Data in R
[2] Raster 02: When Rasters Don't Line Up - Reproject Raster Data in R | NSF NEON | Open Data to Understand our Ecosystems

By Perplexity at For raster object layer in the R programming language, what is the function to set or change crs?

[moniqa]

Thank you for the answer I will try somehow to perform this., maybe with some other Raster Layer, as the layer ttops_chm_p2r is not Raster Layer and that is why I can not use this layer to obtain CRS.

[technocrat]

Or you could just set/change them both to the same crs.

[moniqa]

I did it and it still do not work

chm_p2r_2180 = projectRaster(chm_p2r, crs = 2180)
st_crs(ttops_chm_p2r) = 2180
crowns = dalponte2016(chm_p2r_2180, ttops_chm_p2r)()

Error in geos_op2_geom("intersection", x, y, ...) : 
  st_crs(x) == st_crs(y) is not TRUE

[technocrat]

@moniqa

chm_p2r_2180

Can you provide a link to an object that is the same class as this?

[moniqa]

You mean the link to download the chm_p2r_2180 file?

[technocrat]

Or one like it in an R dataset. I can’t intuit what it is just by a name

[moniqa]

76503_1214667_M-33-35-C-c-1-4-2-1_norm_CHMp2r.tif and 1 other file

That file (chm_p2r_2180) is virtual raster full of TIFs.

I checked both chm_p2r and ttops_chm_p2r in QGIS software - they are in 2180 CRS and on the same place...

[moniqa]

My last post was temporarily hidden, probably because of link with the data.
chm_p2r_2180 is virtual raster. As I have erros with on-disk raster I did chm_p2r = readAll(chm_p2r).
chm_p2r_2180 is full of TIFs.
ttops_chm_p2r is shp.

I checked the CRS on QGIS software and both chm_p2r_2180 and ttops_chm_p2r are in CRS 2180.

[technocrat]

Ping me when it shows up? It's been a long time since I've worked in raster and as a diagnostician I rely heavily on being able to inspect objects.

[technocrat]

@moniqa

ttops_chm_p2r = st_read("ttops.shp")

What is the crs of this object?

library(raster)
#> Loading required package: sp
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)
chm_p2r = raster("/home/roc/projects/demo/rasterize_canopy.vrt")
# it's a RasterLayer class S4 object
str(chm_p2r)
#> Formal class 'RasterLayer' [package "raster"] with 13 slots
#>   ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
#>   .. .. ..@ name        : chr "/home/roc/projects/demo/rasterize_canopy.vrt"
#>   .. .. ..@ datanotation: chr "FLT4S"
#>   .. .. ..@ byteorder   : chr "little"
#>   .. .. ..@ nodatavalue : num -Inf
#>   .. .. ..@ NAchanged   : logi FALSE
#>   .. .. ..@ nbands      : int 1
#>   .. .. ..@ bandorder   : chr "BIL"
#>   .. .. ..@ offset      : int 0
#>   .. .. ..@ toptobottom : logi TRUE
#>   .. .. ..@ blockrows   : Named int 128
#>   .. .. .. ..- attr(*, "names")= chr "rows"
#>   .. .. ..@ blockcols   : Named int 128
#>   .. .. .. ..- attr(*, "names")= chr "cols"
#>   .. .. ..@ driver      : chr "gdal"
#>   .. .. ..@ open        : logi FALSE
#>   ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
#>   .. .. ..@ values    : logi(0) 
#>   .. .. ..@ offset    : num 0
#>   .. .. ..@ gain      : num 1
#>   .. .. ..@ inmemory  : logi FALSE
#>   .. .. ..@ fromdisk  : logi TRUE
#>   .. .. ..@ isfactor  : logi FALSE
#>   .. .. ..@ attributes: list()
#>   .. .. ..@ haveminmax: logi TRUE
#>   .. .. ..@ min       : num -1.81
#>   .. .. ..@ max       : num 39.8
#>   .. .. ..@ band      : int 1
#>   .. .. ..@ unit      : chr ""
#>   .. .. ..@ names     : chr "rasterize_canopy"
#>   ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
#>   .. .. ..@ type      : chr(0) 
#>   .. .. ..@ values    : logi(0) 
#>   .. .. ..@ color     : logi(0) 
#>   .. .. ..@ names     : logi(0) 
#>   .. .. ..@ colortable: logi(0) 
#>   ..@ title   : chr(0) 
#>   ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
#>   .. .. ..@ xmin: num 346458
#>   .. .. ..@ xmax: num 372613
#>   .. .. ..@ ymin: num 354461
#>   .. .. ..@ ymax: num 373968
#>   ..@ rotated : logi FALSE
#>   ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
#>   .. .. ..@ geotrans: num(0) 
#>   .. .. ..@ transfun:function ()  
#>   ..@ ncols   : int 26155
#>   ..@ nrows   : int 19507
#>   ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
#>   .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
#>   .. .. ..$ comment: chr "BOUNDCRS[\n    SOURCECRS[\n        PROJCRS[\"unknown\",\n            BASEGEOGCRS[\"unknown\",\n                "| __truncated__
#>   ..@ srs     : chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
#>   ..@ history : list()
#>   ..@ z       : list()
# it has crs of
chm_p2r@crs@projargs
#> [1] "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

^{Created on 2023-06-28 with reprex v2.0.2}

[moniqa]

According to QGIS software both of them are in CRS 2180, so ttops_chm_p2r = st_read("ttops.shp") is in CRS 2180.

[technocrat]

@moniqa

CRS 2180

That's an abbreviation only.

library(raster)
library(sf)

r <- raster("/home/roc/projects/demo/rasterize_canopy.vrt")
r@crs@projargs

# read in YOUR shape file, however you do it
s <- st_read(system.file("shape/nc.shp", package="sf"))
# use compareCRS() to make sure that both the objects are FULLY identical``
compareCRS(r,s)
# set the crs equal to the crs of r
s <- st_set_crs(s, r)
check again
compareCRS(r,s)

[moniqa]

What does it mean?

That's an abbreviation only.

This is how it looks like:

> compareCRS(chm_p2r,ttops_chm_p2r_wj)
[1] TRUE
> ttops_chm_p2r_wj = st_set_crs(ttops_chm_p2r_wj, chm_p2r)
Error in if (is.na(x)) NA_crs_ else if (inherits(x, "crs")) x else if (is.character(x)) { : 
  argument is not interpretable as logical
> rasterOptions(todisk = TRUE)
> korony_wj = dalponte2016(chm_p2r, ttops_chm_p2r_wj)()
Error in geos_op2_geom("intersection", x, y, ...) : 
  st_crs(x) == st_crs(y) is not TRUE

[technocrat]

A crs looks not like "CRS 2180" but like

+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

@moniqa

ttops_chm_p2r = st_read("ttops.shp")
ttops_chm_p2r = st_cast(ttops_chm_p2r, to = "POINT")

Try it with both versions?

[moniqa]

I did ttops_chm_p2r = st_cast(ttops_chm_p2r, to = "POINT")after ttops_chm_p2r = st_read("ttops.shp"), because when it was only ttops_chm_p2r = st_read("ttops.shp") there was an error. But I will try once again, as I use this script for several data sets and for some of them all the script works and for some do not.

EDIT:
After only ttops_chm_p2r = st_read("ttops.shp") and crowns = dalponte2016(chm_p2r, ttops_chm_p2r_wj)()there is an Error: treetops must be a spatial object of points