Hello, I would like to determines the intersection points between a great circle and a small circle.
In R, I only found gcIntersect function in geosphere package for the intersections of two great circles which is ready to use. Did a lot of search, found gcxsc function in octave forgepackage mapping. https://octave.sourceforge.io/mapping/function/gcxsc.html And use https://github.com/kvasilopoulos/octaver package to call and communication octave in R, but it is slowly. Is there any ready to use functon in R package for determines the intersection points between a great circle and a small circle I missed？ Thanks a lot. Kind regards, Qingshan Zhao qingshanz...@outlook.com From: r-sig-geo-request Date: 2021-10-21 18:00 To: r-sig-geo Subject: R-sig-Geo Digest, Vol 218, Issue 10 Send R-sig-Geo mailing list submissions to r-sig-geo@r-project.org To subscribe or unsubscribe via the World Wide Web, visit https://stat.ethz.ch/mailman/listinfo/r-sig-geo or, via email, send a message with subject or body 'help' to r-sig-geo-requ...@r-project.org You can reach the person managing the list at r-sig-geo-ow...@r-project.org When replying, please edit your Subject line so it is more specific than "Re: Contents of R-sig-Geo digest..." Today's Topics: 1. question on raster Moran's I statistical significance (Gabriel Cotlier) 2. Re: question on raster Moran's I statistical significance (Roger Bivand) 3. Re: question on raster Moran's I statistical significance (Roger Bivand) 4. Re: question on raster Moran's I statistical significance (Gabriel Cotlier) 5. Re: question on raster Moran's I statistical significance (Roger Bivand) ---------------------------------------------------------------------- Message: 1 Date: Wed, 20 Oct 2021 16:20:45 +0300 From: Gabriel Cotlier <gabikl...@gmail.com> To: r-sig-geo <r-sig-geo@r-project.org> Subject: [R-sig-Geo] question on raster Moran's I statistical significance Message-ID: <CAAKwTDHO1rah2NX+z2p83kdypX+XHt=4fYSXT=SJixM=igj...@mail.gmail.com> Content-Type: text/plain; charset="utf-8" Hello, I would like to estimate the Moran's *I* coefficient for raster data and together with the statical significance of the spatial autocorrelation obtained. I found that the raster package function Moran() although calculates the spatial autocorrelation index it apparently does not give directly the statical significance of the results obtained : https://search.r-project.org/CRAN/refmans/raster/html/autocor.html Could it be be possible to obtain the statistical significance of the results with either raster package or similar one? Thanks a lot. Kind regards, Gabriel [[alternative HTML version deleted]] ------------------------------ Message: 2 Date: Wed, 20 Oct 2021 19:42:28 +0200 (CEST) From: Roger Bivand <roger.biv...@nhh.no> To: Gabriel Cotlier <gabikl...@gmail.com> Cc: r-sig-geo <r-sig-geo@r-project.org> Subject: Re: [R-sig-Geo] question on raster Moran's I statistical significance Message-ID: <91dcba60-909f-ed29-1f4-f39c63...@reclus2.nhh.no> Content-Type: text/plain; charset="us-ascii"; Format="flowed" On Wed, 20 Oct 2021, Gabriel Cotlier wrote: > Hello, > > I would like to estimate the Moran's *I* coefficient for raster data and > together with the statical significance of the spatial autocorrelation > obtained. > > I found that the raster package function Moran() although calculates the > spatial autocorrelation index it apparently does not give directly the > statical significance of the results obtained : > https://search.r-project.org/CRAN/refmans/raster/html/autocor.html > > Could it be be possible to obtain the statistical significance of the > results with either raster package or similar one? fortunes::fortune("This is R") library(raster) r <- raster(nrows=10, ncols=10) values(r) <- 1:ncell(r) f <- matrix(c(0,1,0,1,0,1,0,1,0), nrow=3) (rI <- Moran(r, f)) r1 <- r nsim <- 499 res <- numeric(nsim) set.seed(1) for (i in 1:nsim) { values(r1) <- values(r)[sample(prod(dim(r)))] res[i] <- Moran(r1, f) } Hope-type tests date back to Cliff and Ord; they are permutation bootstraps. r_g <- as(r, "SpatialPixelsDataFrame") library(spdep) nb <- poly2nb(as(r_g, "SpatialPolygons"), queen=FALSE) set.seed(1) o <- moran.mc(r_g$layer, nb2listw(nb, style="B"), nsim=nsim, return_boot=TRUE) x_a <- range(c(o$t, o$t0, res, rI)) plot(density(o$t), xlim=x_a) abline(v=o$t0) lines(density(res), lty=2) abline(v=rI, lty=2) It is not immediately obvious from the code of raster::Moran() why it is different, possibly because of padding the edges of the raster and thus increasing the cell count. For added speed, the bootstrap can be parallelized in both cases; polygon boundaries are perhaps not ideal. Hope this clarifies. Always provide a reproducible example, never post HTML mail. Roger Bivand > > Thanks a lot. > Kind regards, > Gabriel > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Roger Bivand Emeritus Professor Department of Economics, Norwegian School of Economics, Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. e-mail: roger.biv...@nhh.no https://orcid.org/0000-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en ------------------------------ Message: 3 Date: Wed, 20 Oct 2021 19:45:01 +0200 (CEST) From: Roger Bivand <roger.biv...@nhh.no> To: Gabriel Cotlier <gabikl...@gmail.com> Cc: r-sig-geo <r-sig-geo@r-project.org> Subject: Re: [R-sig-Geo] question on raster Moran's I statistical significance Message-ID: <3e639398-e137-5546-1848-5a3afb809...@reclus2.nhh.no> Content-Type: text/plain; charset="us-ascii"; Format="flowed" On Wed, 20 Oct 2021, Gabriel Cotlier wrote: > Hello, > > I would like to estimate the Moran's *I* coefficient for raster data and > together with the statical significance of the spatial autocorrelation > obtained. > > I found that the raster package function Moran() although calculates the > spatial autocorrelation index it apparently does not give directly the > statical significance of the results obtained : > https://search.r-project.org/CRAN/refmans/raster/html/autocor.html > > Could it be be possible to obtain the statistical significance of the > results with either raster package or similar one? fortunes::fortune("This is R") library(raster) r <- raster(nrows=10, ncols=10) values(r) <- 1:ncell(r) f <- matrix(c(0,1,0,1,0,1,0,1,0), nrow=3) (rI <- Moran(r, f)) r1 <- r nsim <- 499 res <- numeric(nsim) set.seed(1) for (i in 1:nsim) { values(r1) <- values(r)[sample(prod(dim(r)))] res[i] <- Moran(r1, f) } Hope-type tests date back to Cliff and Ord; they are permutation bootstraps. r_g <- as(r, "SpatialPixelsDataFrame") library(spdep) nb <- poly2nb(as(r_g, "SpatialPolygons"), queen=FALSE) set.seed(1) o <- moran.mc(r_g$layer, nb2listw(nb, style="B"), nsim=nsim, return_boot=TRUE) x_a <- range(c(o$t, o$t0, res, rI)) plot(density(o$t), xlim=x_a) abline(v=o$t0) lines(density(res), lty=2) abline(v=rI, lty=2) It is not immediately obvious from the code of raster::Moran() why it is different, possibly because of padding the edges of the raster and thus increasing the cell count. For added speed, the bootstrap can be parallelized in both cases; polygon boundaries are perhaps not ideal. Hope this clarifies. Always provide a reproducible example, never post HTML mail. Roger Bivand > > Thanks a lot. > Kind regards, > Gabriel > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Roger Bivand Emeritus Professor Department of Economics, Norwegian School of Economics, Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. e-mail: roger.biv...@nhh.no https://orcid.org/0000-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en ------------------------------ Message: 4 Date: Thu, 21 Oct 2021 08:36:51 +0300 From: Gabriel Cotlier <gabikl...@gmail.com> To: roger.biv...@nhh.no Cc: r-sig-geo <r-sig-geo@r-project.org> Subject: Re: [R-sig-Geo] question on raster Moran's I statistical significance Message-ID: <caakwtdhmsathqr3nzlzpwte++t91ood_bvdv8dpg-4uo4hj...@mail.gmail.com> Content-Type: text/plain; charset="utf-8" Hello Thank you very much. I have large raster layers and would like to ask, in order to reduce the processing time of the simulation, choosing a smaller nsim value could help ? and if so, what could be the minimum nsim value recommended? Thanks a lot again. Kind regards, On Wed, Oct 20, 2021 at 8:45 PM Roger Bivand <roger.biv...@nhh.no> wrote: > On Wed, 20 Oct 2021, Gabriel Cotlier wrote: > > > Hello, > > > > I would like to estimate the Moran's *I* coefficient for raster data and > > together with the statical significance of the spatial autocorrelation > > obtained. > > > > I found that the raster package function Moran() although calculates the > > spatial autocorrelation index it apparently does not give directly the > > statical significance of the results obtained : > > https://search.r-project.org/CRAN/refmans/raster/html/autocor.html > > > > Could it be be possible to obtain the statistical significance of the > > results with either raster package or similar one? > > > fortunes::fortune("This is R") > > > library(raster) > r <- raster(nrows=10, ncols=10) > values(r) <- 1:ncell(r) > f <- matrix(c(0,1,0,1,0,1,0,1,0), nrow=3) > (rI <- Moran(r, f)) > r1 <- r > nsim <- 499 > res <- numeric(nsim) > set.seed(1) > for (i in 1:nsim) { > values(r1) <- values(r)[sample(prod(dim(r)))] > res[i] <- Moran(r1, f) > } > > Hope-type tests date back to Cliff and Ord; they are permutation > bootstraps. > > r_g <- as(r, "SpatialPixelsDataFrame") > library(spdep) > nb <- poly2nb(as(r_g, "SpatialPolygons"), queen=FALSE) > set.seed(1) > o <- moran.mc(r_g$layer, nb2listw(nb, style="B"), nsim=nsim, > return_boot=TRUE) > x_a <- range(c(o$t, o$t0, res, rI)) > plot(density(o$t), xlim=x_a) > abline(v=o$t0) > lines(density(res), lty=2) > abline(v=rI, lty=2) > > It is not immediately obvious from the code of raster::Moran() why it is > different, possibly because of padding the edges of the raster and thus > increasing the cell count. > > For added speed, the bootstrap can be parallelized in both cases; polygon > boundaries are perhaps not ideal. > > Hope this clarifies. Always provide a reproducible example, never post > HTML > mail. > > Roger Bivand > > > > > > Thanks a lot. > > Kind regards, > > Gabriel > > > > [[alternative HTML version deleted]] > > > > _______________________________________________ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > > -- > Roger Bivand > Emeritus Professor > Department of Economics, Norwegian School of Economics, > Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. > e-mail: roger.biv...@nhh.no > https://orcid.org/0000-0003-2392-6140 > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en > -- Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] ------------------------------ Message: 5 Date: Thu, 21 Oct 2021 10:40:18 +0200 (CEST) From: Roger Bivand <roger.biv...@nhh.no> To: Gabriel Cotlier <gabikl...@gmail.com> Cc: r-sig-geo <r-sig-geo@r-project.org> Subject: Re: [R-sig-Geo] question on raster Moran's I statistical significance Message-ID: <895fb754-b644-5543-d213-228764873...@reclus2.nhh.no> Content-Type: text/plain; charset="us-ascii"; Format="flowed" On Thu, 21 Oct 2021, Gabriel Cotlier wrote: > Hello > > Thank you very much. > I have large raster layers and would like to ask, in order to reduce the > processing time of the simulation, choosing a smaller nsim value could help > ? and if so, what could be the minimum nsim value recommended? If you examine the code, you see that raster::Moran() always performs multiple raster::focal() calls, effectively constructing the spatial weights on each run. The problem is not nsim, it is the way that raster::Moran() is constructed. Indeed, using a Hope-type test gives the same inferential outcome as using asymptotic methods for global tests for spatial autocorrelation, and is superfluous, but no other approach is feasible for raster::Moran() (which by the way only seems to use 4 neighbours although claiming it uses 8). 1. How large is large? It is possible to generate nb neighbour objects for large data sets, making the use of spdep::moran.test() feasible, certainly for tens of thousands of observations (all the census tracts in coterminous US, for example). This uses distances between raster cell centres to find neighbours: > library(spdep) Loading required package: sp Loading required package: spData Loading required package: sf Linking to GEOS 3.10.0, GDAL 3.3.2, PROJ 8.1.1 > crds <- expand.grid(x=1:800, y=1:1200) > dim(crds) [1] 960000 2 > grd <- st_as_sf(crds, coords=c("x", "y")) > grd$z <- runif(nrow(grd)) > system.time(dnbr <- dnearneigh(grd, 0, 1.01)) user system elapsed 30.065 0.235 30.381 > dnbr Neighbour list object: Number of regions: 960000 Number of nonzero links: 3836000 Percentage nonzero weights: 0.0004162326 Average number of links: 3.995833 > system.time(dnbq <- dnearneigh(grd, 0, 1.42)) user system elapsed 54.502 0.080 54.699 > dnbq Neighbour list object: Number of regions: 960000 Number of nonzero links: 7668004 Percentage nonzero weights: 0.0008320317 Average number of links: 7.987504 Once the neighbour objects are ready, conversion to spatial weights objects takes some time, and computing I with a constant mean model depends on the numbers of neighbours: > system.time(lwr <- nb2listw(dnbr, style="B")) user system elapsed 6.694 0.000 6.722 > system.time(Ir <- moran.test(grd$z, lwr)) user system elapsed 24.371 0.000 24.470 > Ir Moran I test under randomisation data: grd$z weights: lwr Moran I statistic standard deviate = -0.06337, p-value = 0.5253 alternative hypothesis: greater sample estimates: Moran I statistic Expectation Variance -4.679899e-05 -1.041668e-06 5.213749e-07 > system.time(lwq <- nb2listw(dnbq, style="B")) user system elapsed 6.804 0.000 6.828 > system.time(Iq <- moran.test(grd$z, lwq)) user system elapsed 46.703 0.012 46.843 > Iq Moran I test under randomisation data: grd$z weights: lwq Moran I statistic standard deviate = -0.70373, p-value = 0.7592 alternative hypothesis: greater sample estimates: Moran I statistic Expectation Variance -3.604417e-04 -1.041668e-06 2.608222e-07 2. The larger your N, the less likely that the test means anything at all, because the assumption is that the observed entities are not simply arbitrary products of, say, resolution. If you think of global Moran's I as a specification test of a regression of the variable of interest on the constant (the mean model is just the constant), for raster data the resolution controls the outcome (downscaling/upscaling will shift Moran's I). If you include covariates, patterning in the residuals of a richer model may well abate. Hope this clarifies, Roger > Thanks a lot again. > Kind regards, > > > On Wed, Oct 20, 2021 at 8:45 PM Roger Bivand <roger.biv...@nhh.no> wrote: > >> On Wed, 20 Oct 2021, Gabriel Cotlier wrote: >> >>> Hello, >>> >>> I would like to estimate the Moran's *I* coefficient for raster data and >>> together with the statical significance of the spatial autocorrelation >>> obtained. >>> >>> I found that the raster package function Moran() although calculates the >>> spatial autocorrelation index it apparently does not give directly the >>> statical significance of the results obtained : >>> https://search.r-project.org/CRAN/refmans/raster/html/autocor.html >>> >>> Could it be be possible to obtain the statistical significance of the >>> results with either raster package or similar one? >> >> >> fortunes::fortune("This is R") >> >> >> library(raster) >> r <- raster(nrows=10, ncols=10) >> values(r) <- 1:ncell(r) >> f <- matrix(c(0,1,0,1,0,1,0,1,0), nrow=3) >> (rI <- Moran(r, f)) >> r1 <- r >> nsim <- 499 >> res <- numeric(nsim) >> set.seed(1) >> for (i in 1:nsim) { >> values(r1) <- values(r)[sample(prod(dim(r)))] >> res[i] <- Moran(r1, f) >> } >> >> Hope-type tests date back to Cliff and Ord; they are permutation >> bootstraps. >> >> r_g <- as(r, "SpatialPixelsDataFrame") >> library(spdep) >> nb <- poly2nb(as(r_g, "SpatialPolygons"), queen=FALSE) >> set.seed(1) >> o <- moran.mc(r_g$layer, nb2listw(nb, style="B"), nsim=nsim, >> return_boot=TRUE) >> x_a <- range(c(o$t, o$t0, res, rI)) >> plot(density(o$t), xlim=x_a) >> abline(v=o$t0) >> lines(density(res), lty=2) >> abline(v=rI, lty=2) >> >> It is not immediately obvious from the code of raster::Moran() why it is >> different, possibly because of padding the edges of the raster and thus >> increasing the cell count. >> >> For added speed, the bootstrap can be parallelized in both cases; polygon >> boundaries are perhaps not ideal. >> >> Hope this clarifies. Always provide a reproducible example, never post >> HTML >> mail. >> >> Roger Bivand >> >> >>> >>> Thanks a lot. >>> Kind regards, >>> Gabriel >>> >>> [[alternative HTML version deleted]] >>> >>> _______________________________________________ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> -- >> Roger Bivand >> Emeritus Professor >> Department of Economics, Norwegian School of Economics, >> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. >> e-mail: roger.biv...@nhh.no >> https://orcid.org/0000-0003-2392-6140 >> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en >> > > > -- Roger Bivand Emeritus Professor Department of Economics, Norwegian School of Economics, Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway. e-mail: roger.biv...@nhh.no https://orcid.org/0000-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en ------------------------------ Subject: Digest Footer _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo ------------------------------ End of R-sig-Geo Digest, Vol 218, Issue 10 ****************************************** [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
