Should coercion from ppp to sf respect window scaling?

[rsbivand]

For:

library(spatstat)
data(japanesepines)
japanesepines
#Planar point pattern: 65 points
#window: rectangle = [0, 1] x [0, 1] units (one unit = 5.7 metres)
str(japanesepines)
#List of 5
# $ window    :List of 4
# ..$ type  : chr "rectangle"
#  ..$ xrange: num [1:2] 0 1
#  ..$ yrange: num [1:2] 0 1
#  ..$ units :List of 3
#  .. ..$ singular  : chr "metre"
#  .. ..$ plural    : chr "metres"
#  .. ..$ multiplier: num 5.7
#  .. ..- attr(*, "class")= chr "unitname"
#  ..- attr(*, "class")= chr "owin"
# $ n         : int 65
# $ x         : num [1:65] 0.09 0.29 0.38 0.39 0.48 0.59 0.65 0.67 0.73 0.79 ...
# $ y         : num [1:65] 0.09 0.02 0.03 0.18 0.03 0.02 0.16 0.13 0.13 0.03 ...
# $ markformat: chr "none"
# - attr(*, "class")= chr "ppp"

maptools has:

library(maptools)
spjpines0 <- as(japanesepines, "SpatialPoints")
library(sp)
bbox(spjpines0)
#     min max
#[1,]   0 5.7
#[2,]   0 5.7
head(coordinates(spjpines0))
#        mx    my
#[1,] 0.513 0.513
#[2,] 1.653 0.114
#[3,] 2.166 0.171
#[4,] 2.223 1.026
#[5,] 2.736 0.171
#[6,] 3.363 0.114

that is both the bounding box set from the owin object and the coordinates themselves are rescaled by window$units$multiplier. In sf:

library(sf)
oo <- st_as_sf(japanesepines)
oo
#Simple feature collection with 66 features and 1 field
#Geometry type: GEOMETRY
#Dimension:     XY
#Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
#CRS:           NA
#First 10 features:
#    label                           geom
#1  window POLYGON ((0 0, 1 0, 1 1, 0 ...
#2   point              POINT (0.09 0.09)
#3   point              POINT (0.29 0.02)
#4   point              POINT (0.38 0.03)
#5   point              POINT (0.39 0.18)
#6   point              POINT (0.48 0.03)

Found while checking ASDAR ch. 7 for retiring maptools. Was there a reason for not re-scaling in the sf implementation? @rubak is re-scaling actively used with the owin multiplier component? Do we need to ensure round-tripping (maptools/sp do not roundtrip the multiplier)?

[rubak]

This is a bit of a though one. As I see it the multiplier only exists because of historic datasets that have been rescaled to the unit square. In practical matters this should not be very relevant since I don't see any use for it when you read in data and analyze it, but it should of course be handled correctly. It definitely makes sense to do what maptools does, but on the other hand it may be a bit confusing that the coordinates change just because you change from spatstat to sp representation. The spatstat philosophy has been to leave anything related to units to be handled by the user. So maybe if you want the converted units to be in meters you should do the conversion yourself (just type jp <- rescale(japanesepines) and then convert jp).

In conclusion: I think it is better to leave the conversion as is is sf. I will also try to ask @baddstats for his opinion on this.

[baddstats]

I agree with @rubak overall: it's probably better to leave the conversion as it is in sf.

The unit multiplier was introduced in spatstat because we wanted to handle spatial scale, while many of the historic (1970's) datasets had been rescaled to the unit square. But they were rescaled to the unit square for a reason: to avoid numerical problems (overflow, instability, failure to converge). Many of the numerical problems have been overcome in the last 40 years, but there are still many remaining --- cases where the result of analysing a dataset depends on the scaling of the coordinate system. After all, many of the basic tools such as optimisation algorithms (optim, etc) do not automatically compensate for scale. Consequently the policy in spatstat remains that the units will not be changed except by explicit command of the user.

The question is whether this spatstat policy should still be followed in other packages when they convert spatstat data. I would suggest doing this, for the sake of the beginner. To help the beginner, one could issue a message like "unit multiplier discarded"