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

[Michael Sumner]

On Sat, Jul 31, 2021 at 6:09 AM Alexandre Santos <
alexandresantosbr using 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 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
> >
> >
>


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

	[[alternative HTML version deleted]]