Re: [R-sig-Geo] Create a ppp object in UTM with mixture of zones
On Sat, Jul 31, 2021 at 6:09 AM Alexandre Santos <
alexandresanto...@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 <
mdsum...@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@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
Re: [R-sig-Geo] Create a ppp object in UTM with mixture of zones
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
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,
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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Create a ppp object in UTM with mixture of zones
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@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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
