Create a ppp object in UTM with mixture of zones

On Sat, Jul 31, 2021 at 6:09 AM Alexandre Santos <

Thanks, Mike

Problem solved with a transformation of to Lambert azimuthal equal-area

projection centred on longitude and latitude of 0:

# Packages
library(spatstat)
library(sf)
library(sp)
library(raster)


# Open spatial data set in GitHub
sp_ds<-read.csv("

https://raw.githubusercontent.com/Leprechault/trash/master/myspds.csv")

str(sp_ds)
#'data.frame':  4458 obs. of  2 variables:
# $ Lat : num  9.17 9.71 9.12 9.12 9.71 ...
# $ Long: num  35.8 35.5 35.8 35.8 35.5 ...


sp_ds <- st_read("

https://raw.githubusercontent.com/Leprechault/trash/master/myspds.csv")

sp_ds_sf <- st_as_sf(sp_ds, coords = c("Lat", "Long"), crs = 4326)


# Here the Mike tip!!
sp_ds_laea = st_transform(sp_ds_sf,
                           crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=0

+lat_0=0")


  great, but you should centre the projection for your data - +lon_0 and
+lat_0 specify a central-ish point for your longitude,latitude values. Plot
the data on a map in long lat, and then plot them once projected - it looks
pretty weird with centre lon,lat at 0,0 . - you want more like "+proj=laea
+lon_0=29.9 +lat_0=13.4 +datum=WGS84" - the x_0 and y_0 may be specified,
they just move the 0,0 of the projection away from the centre of your map
(nice to not have negative values). It's likely you could find a local
authority projection that is defined for general use (laea, lcc, aeqd are
common).

Also, you are using Lat,Long order which is technically correct for
EPSG:4326, other variants of longitude,latitude specifications use the more
straightforward coordinate order, it's configurable and sf handles it etc.
but you should definitely check :)

Finally, it's probably worth saying in this case for these data you could
probably choose the central (of the 3) UTM zones involved and get just as
good properties ( I never see that said in the wild). But, it's not really
possible to have overarching advice here because your data and your context
and your distances can be a very broad set of cases, that's why different
families of projections exist, for options. You can calculate distance and
area in long,lat with geodist, geosphere, sf and many other packages - and
sf is incorporating s2 to enable that always for data in long,lat - but I
assume you need an actual coordinate system for spatstat to operate in its
own context (so your approach is exactly right). The choice of projection
can be subtle and require expertise, but a local central laea for smallish
regions is pretty fail safe.  Over a small area like this, most projections
will be exactly the same in terms of distance, area, shape  - but not
knowing the broader case of data you might want to deal with means such
advice needs to be issued with care. :)

HTH

# Create boudaries using sf
ch_ds <- st_convex_hull(st_union(sp_ds_laea))
W<- as.owin(ch_ds)


# Create a ppp Point Pattern with sf object
out.ppp <- as.ppp(X=st_coordinates(sp_ds_laea),W=W)
plot(out.ppp)




# Make a CRS test
f1 <- ppm(out.ppp~1)
E <-envelope(f1, Kinhom, nsim = 19, global = TRUE, correction = "border")
plot(E)

Best wishes,

Alexandre








Em sexta-feira, 30 de julho de 2021 15:41:08 BRT, Michael Sumner <

mdsumner at gmail.com> escreveu:

Don't use UTM, it's a legacy of a bygone era. Use a local custom

projection such as laea:

https://twitter.com/mdsumner/status/1136794870113218561?s=19

I know it's common advice, utm for projection choice but it's trash

propagated by folks who should know better. Traversing zones for projection
choice for basic distance properties is unnecessary.

Best, Mike

On Fri, 30 Jul 2021, 08:21 Alexandre Santos via R-sig-Geo, <

Dear Members,


I'd like to use my spatial data set ("sp_ds") as ppp (spatstat Point

Pattern object), an I try to do:

# Open spatial data set in GitHub
library(spatstat)
library(sf)
library(sp)


sp_ds<-read.csv("

https://raw.githubusercontent.com/Leprechault/trash/master/myspds.csv")

str(sp_ds)
#'data.frame':  4458 obs. of  2 variables:
# $ Lat : num  9.17 9.71 9.12 9.12 9.71 ...
# $ Long: num  35.8 35.5 35.8 35.8 35.5 ...


# Create boudaries using sf
sfds = st_as_sf(sp_ds, coords=c("Long","Lat"), crs=4326)
traps<-sp_ds
ch <- chull(traps[,c(2,1)])
coords <- traps[c(ch, ch[1]), ]
coordinates(coords) <- c("Long","Lat")
proj4string(coords) <- CRS("+init=epsg:4326")
W <- owin(poly=cbind(coordinates(coords)[,2],coordinates(coords)[,1]))



# Create a ppp Point Pattern object
out.ppp<-ppp(x=sp_ds$Lat,y=sp_ds$Long,window=W)
plot(out.ppp)


f1 <- ppm(out.ppp~1)
E <-envelope(f1, Kinhom, nsim = 19, global = TRUE, correction =

"border")

plot(E)


#But I'd like to r distance (x axis) in kilometers and for this I need

to convert the coordinate reference system of 4326 to

UTM, a have 3 UTM zones:
#34N bounds: (18.0, 0.0, 24.0, 84.0)
#35N bounds: (24.0, 0.0, 30.0, 84.0)
#36N bounds: (30.0, 0.0, 36.0, 84.0)


Please the area a simple way to create a ppp object in UTM with a

mixture of zones?

Thanks in advance,


Alexandre

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsumner at gmail.com