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
Create a ppp object in UTM with mixture of zones
4 messages · ASANTOS, Michael Sumner
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 at 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 at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
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 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, <r-sig-geo at 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 at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
On Sat, Jul 31, 2021 at 6:09 AM Alexandre Santos <
alexandresantosbr at yahoo.com.br> wrote:
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, <
r-sig-geo at 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 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