[R-sig-Geo] Create a ppp object in UTM with mixture of zones

[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")
# 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 using 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, <r-sig-geo using r-project.org> wrote:
> 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 using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
>