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

moniqa · June 27, 2023, 8:29pm

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 · June 28, 2023, 1:56am

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:

Load the raster package:

library(raster)

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.
Use the crs() function to obtain the CRS of raster1:

raster1_crs <- crs(raster1)

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 · June 28, 2023, 8:06am

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 · June 28, 2023, 8:08am

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

moniqa · June 28, 2023, 8:26am

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 · June 28, 2023, 8:29am

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

moniqa · June 28, 2023, 8:32am

You mean the link to download the chm_p2r_2180 file?

technocrat · June 28, 2023, 8:37am

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

moniqa · June 28, 2023, 8:46am

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 · June 28, 2023, 9:57am

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 · June 28, 2023, 10:14am

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 · June 29, 2023, 6:10am

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 · June 29, 2023, 6:47am

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

technocrat · June 29, 2023, 9:01am

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 · June 29, 2023, 4:46pm

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 · June 29, 2023, 11:45pm

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

Try it with both versions?

moniqa · June 30, 2023, 6:56am

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

system · July 21, 2023, 6:01am

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.