Set CRS for an arbitrary Euclidean plane

[aloboa]

Is it possible to set an arbitrary Euclidean plane as CRS in terra?
In some cases (eg. relatively small field plots) there is really no standard CRS.
According to
https://gis.stackexchange.com/questions/27699/how-to-represent-an-imaginary-flat-world-in-qgis-without-using-a-round-world-crs
I could use eg. World Mercator, but QGIS has the option of defining No projection by default.
At the moment, any raster with no CRS information is assumed to be lon/lat WGS 84 in terra.

[brownag]

I believe you can set crs="" argument, or use crs(x) <- "" (with custom x/y min/max arguments to set extent as needed) if you do not want the defaults that correspond to WGS84

[kadyb]

I assume @aloboa is referring to raster data, but for vector data it would probably be better to set up some CRS because some operations may not work.

library("terra")
v = vect(cbind(0, 0), type = "points", crs = "")
buffer(v, 1)
#> Error: [buffer] crs not defined

[brownag]

Could use a LOCAL_CS or ENGCRS WKT string to get more specific and reoslve case brought up by @kadyb. Something like:

vect(cbind(0, 0), type = "points", crs = 'LOCAL_CS["Cartesian (Meter)",LOCAL_DATUM["Local Datum",0],UNIT["Meter",1.0],AXIS["X",NORTH],AXIS["Y",EAST]]')

[aloboa]

Could use a LOCAL_CS or ENGCRS WKT string to get more specific and reoslve case brought up by @kadyb. Something like:

vect(cbind(0, 0), type = "points", crs = 'LOCAL_CS["Cartesian (Meter)",LOCAL_DATUM["Local Datum",0],UNIT["Meter",1.0],AXIS["X",NORTH],AXIS["Y",EAST]]')

This is actually better, as the "no projection" alternative does not always work (even in QGIS, you can set "No projection (or Unknown/non-Earth projection" in Options/CRS but cannot, eg, create an object with no projection).

[aloboa]

@brownag, why "AXIS["X",NORTH],AXIS["Y",EAST]]", should it not be AXIS["Y",NORTH],AXIS["X",EAST]] ?

[brownag]

@brownag, why "AXIS["X",NORTH],AXIS["Y",EAST]]", should it not be AXIS["Y",NORTH],AXIS["X",EAST]] ?

Sorry, that's a typo! No reason

[rhijmans]

Thanks for the useful ideas. I have made a few changes to make this easier to use. You can set the crs that @brownag suggests by using the value "local" (see below). I have also made some changes that avoid errors when this crs is used; but please let me know if you find more cases where this can be done.

SpatRaster

library(terra)
r <- rast(crs="local")
r
#class       : SpatRaster 
#dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter) 

For area computations you need to use "transform=FALSE"

expanse(r, transform=FALSE)
#[1] 64800

cellSize(r, transform=F)
#class       : SpatRaster 
#dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter) 
#source      : memory 
#name        : area 
#min value   :    1 
#max value   :    1 

SpatVector

v <- vect(cbind(c(1,1,2,2,1), c(1,2,2,1,1)), "polygon", crs="local") 
v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 0  (geometries, attributes)
# extent      : 1, 2, 1, 2  (xmin, xmax, ymin, ymax)
# coord. ref. : Cartesian (Meter) 


expanse(v, transform=F)
#[1] 1

[aloboa]

Thanks for taking this suggestion into consideration.
When is this improvement supposed to be included in the regular CRAN version?

[rhijmans]

I expect this to get to CRAN sometime this month. But you can use
install.packages('terra', repos='https://rspatial.r-universe.dev') is installation on Windows or Mac is an issue.

[aloboa]

Given

library(terra)
r <- rast(nrows=10, ncols=5,nlyrs=1, res=1,crs="local")

I get:

> r
class       : SpatRaster 
dimensions  : 180, 360, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : Cartesian (Meter) 
> dim(r)
[1] 180 360   1

Is this correct? I expected the dimension to be 10,5,1 as in

> r <- rast(nrows=10, ncols=5,nlyrs=1, res=1,crs="local",extent=c(0,5,0,10))
> r
class       : SpatRaster 
dimensions  : 10, 5, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 5, 0, 10  (xmin, xmax, ymin, ymax)
coord. ref. : Cartesian (Meter) 
> dim(r)
[1] 10  5  1

[rhijmans]

Perhaps this should return an error:

r <- rast(nrows=10, ncols=5,nlyrs=1, res=1,crs="local")

Because, given the default extent, the combination of rows/columns and resolution is not possible.

[rhijmans]

The manual says:

resolution: numeric vector of length 1 or 2 to set the resolution (see res). If this argument is used, arguments ncol and nrow are ignored

[aloboa]

I wrongly assumed that in this case

r <- rast(nrows=10, ncols=5,nlyrs=1, res=1,crs="local")

the extent was automatically set to extent=c(0,5,0,10). But I understand you must have good reasons for
not doing it. In that case, I agree the error you mention should be returned. Also, an example with the full syntax
to create a raster for an arbitrary field plot would be very appreciated in the documentation, eg
(please correct me if I'm wrong):

to define a plot of 5m x 10m and 1m resolution:

> r <- rast(nrows=10, ncols=5,nlyrs=1, res=1,crs="local",extent=c(0,5,0,10))
> r
class       : SpatRaster 
dimensions  : 10, 5, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 5, 0, 10  (xmin, xmax, ymin, ymax)
coord. ref. : Cartesian (Meter) 

and
for a plot of 5m x 10m and 1 cm resolution

> r <- rast(nrows=1000, ncols=500,nlyrs=1, res=0.01,crs="local",extent=c(0,5,0,10))
> r
class       : SpatRaster 
dimensions  : 1000, 500, 1  (nrow, ncol, nlyr)
resolution  : 0.01, 0.01  (x, y)
extent      : 0, 5, 0, 10  (xmin, xmax, ymin, ymax)
coord. ref. : Cartesian (Meter) 

[rhijmans]

you must have good reasons for not doing it

There are default values for xmin, xmax, ymin, ymax, so it would be unexpected to change these based on the value of other arguments.

Your use of rast has some redundancy. You can specify nrows/ncols or the resolution. If you do both, the nrows/ncols arguments are ignored anyway. There is a slight difference between these two options. If you use resolution, the extent may be adjusted if need be, in order to get an integer number of rows and columns.

For example:

rast(res=1.75, crs="local", extent=c(0,5,0,10))
#class       : SpatRaster 
#dimensions  : 6, 3, 1  (nrow, ncol, nlyr)
#resolution  : 1.75, 1.75  (x, y)
#extent      : 0, 5.25, 0, 10.5  (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter) 

[aloboa]

see
dieghernan/tidyterra#64