Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 1 hour 45 min ago

Re: Intensity Estimation Methods

Sun, 12/17/2017 - 04:44
Adrian and Rolf.
Adrian's second example will solve my problem I guess.
Because I do not know the original intensity of the pattern. I only have an observed point pattern.

Thanks a lot.
Cenk


From: Adrian Baddeley [mailto:[hidden email]]
Sent: Saturday, December 16, 2017 11:58 AM
To: Cenk ���Z <[hidden email]>; Rolf Turner <[hidden email]>
Cc: [hidden email]; Ege Rubak <[hidden email]>
Subject: Re: [R-sig-Geo] Intensity Estimation Methods




Cenk ���Z <[hidden email]<mailto:[hidden email]>> writes:



 I have a spatial point pattern . I am trying to estimate its intensity
> both with a fixed bandwidth and  with an adaptive bandwidth. How could
> I compare the goodness of these two fits? I mean are there any things
> like mse, aic or any other criteria??? I want to compare the
> difference between the estimated intensity and the original pattern's
> intensity.



If you know the 'true' intensity then you could compute, for example, the integrated squared error.

In the spatstat package, if 'lamtrue' is the true intensity and 'lamest' the estimated intensity, given as pixel images (class 'im') then you can just do

     ISE <- integral((lamtrue-lamest)^2)



Alternatives include the Kulback-Leibler divergence

     KL <- integral(log(lamest/lamtrue) * lamtrue)

and the total variation distance

     TV <- integral(abs(lamtrue-lamest))/2



However if you have two competing estimates of the intensity of an observed point pattern, you're probably best to use the point process likelihood. Suppose lam1 and lam2 are pixel images giving  two competing estimates of the intensity of the same point pattern X. Then you could do



      lik1 <- sum(log(lam1[X])) - integral(lam1)

      lik2 <- sum(log(lam2[X])) - integral(lam2)



and compare the likelihoods.



Adrian Baddeley



________________________________
From: Cenk ���Z <[hidden email]<mailto:[hidden email]>>
Sent: Friday, 15 December 2017 8:23 PM
To: Rolf Turner
Cc: [hidden email]<mailto:[hidden email]>; Adrian Baddeley; Ege Rubak
Subject: RE: [R-sig-Geo] Intensity Estimation Methods

Thanks a lot.
I did superimpose the original pattern to pixel images of intensities. With adaptive smothing intensity higher zones are too narrowed.
Also the differences of intensities are getting higher in the study region. The fixed bandwidth choosen with bw.ppl( ) in spatstat give me a better Picture. This is my opinion and also this is in my case . I use this bandwidth as a global bandwidth for adaptive smoothing.


-----Original Message-----
From: Rolf Turner [mailto:[hidden email]]
Sent: Thursday, December 14, 2017 11:14 PM
To: Cenk ���Z <[hidden email]<mailto:[hidden email]>>
Cc: [hidden email]<mailto:[hidden email]>; [hidden email]<mailto:[hidden email]>; Ege Rubak <[hidden email]<mailto:[hidden email]>>
Subject: Re: [R-sig-Geo] Intensity Estimation Methods

On 15/12/17 01:01, Cenk ���Z via R-sig-Geo wrote:
> Hello list,,
>
> I have a spatial point pattern . I am trying to estimate its intensity
> both with a fixed bandwidth and  with an adaptive bandwidth. How could
> I compare the goodness of these two fits? I mean are there any things
> like mse, aic or any other criteria??? I want to compare the
> difference between the estimated intensity and the original pattern's
> intensity.

I think that the following fortune (fortunes::fortune(340)) might be
relevant:

> Bandwidth selection is an unresolved (and possibly unsolvable) problem
> in smoothing, so you're perfectly justified in trying/choosing an
> arbitrary value if it produces good pictures!
>    -- Adrian Baddeley (answering a user's question about the choice of
>       smoothing parameter when using the density.ppp() function from the
>       spatstat package)
>       private communication (March 2013)

I am cc-ing to the man himself to see if he has further comment.

cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276


________________________________

Bu elektronik posta ve onunla iletilen b�t�n dosyalar sadece yukar�da isimleri belirtilen ki�iler aras�nda �zel haberle�me amac�n� ta��makta olup g�nderici taraf�ndan al�nmas� ama�lanan yetkili ger�ek ya da t�zel ki�inin kullan�m�na aittir. E�er bu elektronik posta size yanl��l�kla ula�m��sa, elektronik postan�n i�eri�ini a��klaman�z, kopyalaman�z, y�nlendirmeniz ve kullanman�z kesinlikle yasakt�r. Bu durumda, l�tfen mesaj� geri g�nderiniz ve sisteminizden siliniz. Anadolu �niversitesi bu mesaj�n i�erdi�i bilgilerin do�rulu�u veya eksiksiz oldu�u konusunda herhangi bir garanti vermemektedir. Bu nedenle bu bilgilerin ne �ekilde olursa olsun i�eri�inden, iletilmesinden, al�nmas�ndan ve saklanmas�ndan sorumlu de�ildir. Bu mesajdaki g�r��ler yaln�zca g�nderen ki�iye aittir ve Anadolu �niversitesinin g�r��lerini yans�tmayabilir.

This electronic mail and any files transmitted with it are intended for the private use of the people named above. If you are not the intended recipient and received this message in error, forwarding, copying or use of any of the information is strictly prohibited. Any dissemination or use of this information by a person other than the intended recipient is unauthorized and may be illegal. In this case, please immediately notify the sender and delete it from your system. Anadolu University does not guarantee the accuracy or completeness of any information included in this message. Therefore, by any means Anadolu University is not responsible for the content of the message, and the transmission, reception, storage, and use of the information. The opinions expressed in this message only belong to the sender of it and may not reflect the opinions of Anadolu University.

        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Sat, 12/16/2017 - 03:10
Christopher W. Ryan writes:


> I have a SpatialPolygonsDataFrame of census tract poverty levels in 3
> contiguous counties in the US, called sremsPoverty. I want to use this
> as a predictor in a ppm model. The window for the point pattern is the
> three counties--so an irregular polygon--called sremsWindow.
>
> I understand ppm predictors need to be an image, a tesselation, a funxy,
> a window, or a single number. So I proceed as follows:
>
> ### Poverty
> p <- slot(sremsPoverty, "polygons")
> v <- lapply(p, function(z) { SpatialPolygons(list(z)) })
> sat <- tess(tiles = lapply(v, as.owin) )
> pov.f <- as.function.tess(sat, values = sremsPoverty@data$propPoverty)
>
> Thus pov.f is a spatial function (funxy) that I can use in a call to ppm():
>
> m1 <- ppm(unique.cases.ppp ~ pov.f)
>
> pov.f looks as expected when I plot it. But examing the structure of
> as.im(pov.f) it appears it is a 128 x 128 pixel array, with the value of
> the function at all pixels outside of the irregular polygonal window,
> but inside the bounding rectangle, set to NA.
> I think this is the cause of the NA values I am seeing among the residuals from m1,
> and those NA residuals are causing me some trouble with model diagnostics such as
> rhohat().
> How do I constrain the funxy (or the image I can derive from it) to the
> irregular polygonal window, so as to eliminate the NA values outside the
> window but inside the bounding rectangle? Or can I constrain the
> modeling activity of ppm() to the window?

When you convert the tessellation 'sat' into the function 'pov.f' using as.function.tess, the domain of the function is the union of all the tiles in the tessellation 'sat'. After all, this function pov.f(x) works by figuring out which tile of the tessellation the given location x falls inside, and returns the poverty value associated with that tile. So this function is restricted to the union of the tiles given.

When you fit a model to the point pattern 'unique.cases.ppp', the fitting procedure has to place sample points all over the window associated with the point pattern, and evaluate the covariate 'pov.f' at all these sample points. So this window should not be larger than the window where the poverty values are defined - if it is, then you'll get NA's.

I suggest you do something like
        X <- unique.cases.ppp
        W <- intersect.owin(Window(X), Window(pov.f))
        Y <- X[W]
        m1 <- ppm(Y ~ pov.f)

Adrian Baddeley





        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Intensity Estimation Methods

Sat, 12/16/2017 - 02:58

Cenk ���Z <[hidden email]> writes:


 I have a spatial point pattern . I am trying to estimate its intensity
> both with a fixed bandwidth and  with an adaptive bandwidth. How could
> I compare the goodness of these two fits? I mean are there any things
> like mse, aic or any other criteria??? I want to compare the
> difference between the estimated intensity and the original pattern's
> intensity.


If you know the 'true' intensity then you could compute, for example, the integrated squared error.

In the spatstat package, if 'lamtrue' is the true intensity and 'lamest' the estimated intensity, given as pixel images (class 'im') then you can just do

     ISE <- integral((lamtrue-lamest)^2)


Alternatives include the Kulback-Leibler divergence

     KL <- integral(log(lamest/lamtrue) * lamtrue)

and the total variation distance

     TV <- integral(abs(lamtrue-lamest))/2


However if you have two competing estimates of the intensity of an observed point pattern, you're probably best to use the point process likelihood. Suppose lam1 and lam2 are pixel images giving  two competing estimates of the intensity of the same point pattern X. Then you could do


      lik1 <- sum(log(lam1[X])) - integral(lam1)

      lik2 <- sum(log(lam2[X])) - integral(lam2)


and compare the likelihoods.


Adrian Baddeley


________________________________
From: Cenk ���Z <[hidden email]>
Sent: Friday, 15 December 2017 8:23 PM
To: Rolf Turner
Cc: [hidden email]; Adrian Baddeley; Ege Rubak
Subject: RE: [R-sig-Geo] Intensity Estimation Methods

Thanks a lot.
I did superimpose the original pattern to pixel images of intensities. With adaptive smothing intensity higher zones are too narrowed.
Also the differences of intensities are getting higher in the study region. The fixed bandwidth choosen with bw.ppl( ) in spatstat give me a better Picture. This is my opinion and also this is in my case . I use this bandwidth as a global bandwidth for adaptive smoothing.


-----Original Message-----
From: Rolf Turner [mailto:[hidden email]]
Sent: Thursday, December 14, 2017 11:14 PM
To: Cenk ���Z <[hidden email]>
Cc: [hidden email]; [hidden email]; Ege Rubak <[hidden email]>
Subject: Re: [R-sig-Geo] Intensity Estimation Methods

On 15/12/17 01:01, Cenk ���Z via R-sig-Geo wrote:
> Hello list,,
>
> I have a spatial point pattern . I am trying to estimate its intensity
> both with a fixed bandwidth and  with an adaptive bandwidth. How could
> I compare the goodness of these two fits? I mean are there any things
> like mse, aic or any other criteria??? I want to compare the
> difference between the estimated intensity and the original pattern's
> intensity.

I think that the following fortune (fortunes::fortune(340)) might be
relevant:

> Bandwidth selection is an unresolved (and possibly unsolvable) problem
> in smoothing, so you're perfectly justified in trying/choosing an
> arbitrary value if it produces good pictures!
>    -- Adrian Baddeley (answering a user's question about the choice of
>       smoothing parameter when using the density.ppp() function from the
>       spatstat package)
>       private communication (March 2013)

I am cc-ing to the man himself to see if he has further comment.

cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276


________________________________

Bu elektronik posta ve onunla iletilen b�t�n dosyalar sadece yukar�da isimleri belirtilen ki�iler aras�nda �zel haberle�me amac�n� ta��makta olup g�nderici taraf�ndan al�nmas� ama�lanan yetkili ger�ek ya da t�zel ki�inin kullan�m�na aittir. E�er bu elektronik posta size yanl��l�kla ula�m��sa, elektronik postan�n i�eri�ini a�klaman�z, kopyalaman�z, y�nlendirmeniz ve kullanman�z kesinlikle yasakt�r. Bu durumda, l�tfen mesaj� geri g�nderiniz ve sisteminizden siliniz. Anadolu �niversitesi bu mesaj�n i�erdi�i bilgilerin do�rulu�u veya eksiksiz oldu�u konusunda herhangi bir garanti vermemektedir. Bu nedenle bu bilgilerin ne �ekilde olursa olsun i�eri�inden, iletilmesinden, al�nmas�ndan ve saklanmas�ndan sorumlu de�ildir. Bu mesajdaki g�r��ler yaln�zca g�nderen ki�iye aittir ve Anadolu �niversitesinin g�r��lerini yans�tmayabilir.

This electronic mail and any files transmitted with it are intended for the private use of the people named above. If you are not the intended recipient and received this message in error, forwarding, copying or use of any of the information is strictly prohibited. Any dissemination or use of this information by a person other than the intended recipient is unauthorized and may be illegal. In this case, please immediately notify the sender and delete it from your system. Anadolu University does not guarantee the accuracy or completeness of any information included in this message. Therefore, by any means Anadolu University is not responsible for the content of the message, and the transmission, reception, storage, and use of the information. The opinions expressed in this message only belong to the sender of it and may not reflect the opinions of Anadolu University.

        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Fri, 12/15/2017 - 23:13
On 16/12/17 11:43, Christopher W. Ryan wrote:
> Thanks Rolf. I'm going to have to reflect more on my code and my data,
> to understand better what is going on.
>
> Obviously this won't help you much, without having access to all my data
> and preceeding code,

Too true!

> but the error message that is tripping me up is:
>
>> rhohat(m12, pov.f)
> Error:  the fitted intensity contains NA values

No idea what is causing that to happen, and it's impossible to work it
out without having your data.  (Point pattern and the predictor image.)

I can't find that error message anywhere in the spatstat code.  Is it
possible that you have some other rhohat() function (maybe from some
other package) hanging around and getting in the way?  Using traceback()
*might* provide some insight.

>
> And yet,
>
> table(is.na(fitted(m12)))
> FALSE
>    876
>
> The predicted intensity, however, contains many NA values:
>
> table(is.na(predict(m12)$v))
> FALSE  TRUE
> 10379  6005
This is a *COMPLETE RED HERRING* !!!

Of course you will get NAs from this.  The predicted intensity is an
*image* on the window of the point pattern to which the model is fitted.
All pixels within the bounding box of that window, but outside of the
actual window, have pixel value NA.
>
> I try to force predictions only within my window by specifying locations
> (which I think requires a binary mask window) but get the same result:
>
>> table(is.na(predict(m12, locations = as.mask(sremsWindow))$v))
> FALSE  TRUE
> 10379  6005

No, no, no!!!  The predictions are always "within your window".  That's
*why* they are NA at pixels outside that window.

>
> Does rhohat use fitted values (at quadrature points) or predicted values
> (on a 128 x 128 pixel grid within the window)? Top of p. 415 in your
> book Spatial Point Patterns: Methodology and Applications wtih R seems
> to indicate the latter, while the error message from my rhohat command
> above speaks of fitted values.

I don't understand where that error message is coming from, so I don't
get what it is on about, but essentially rhohat() looks at values of the
predictive covariate and of the fitted intensity at all pixel centres
within the window.  The 128 x 128 grid is the default here, but can be
changed (in the call to rhohat()).

> And how is a rectangular 128 x 128 grid
> fit in an irregularly-shaped polygonal window?  Maybe that's how NA
> predicted values arise--pixels outside an intra-window rectangular grid
> but still inside the window?

The rectangular grid is over the bounding box of the window.  (Which is
rectangular!!!) Only pixels whose centres are within the window have
non-NA values.

>
>
> And I can see now that no residuals from the model are missing:
>
>> table(is.na(residuals(m12)$val))
> FALSE
>    876
>
> All the NA's in my predicted values *around* my window, but inside the
> bounding rectangle, led me down the garden path.
Indeed.

>
> The origin of most of my predictors, such as pov.f, are shapefiles from
> the US Census Bureau, with a discrete value of poverty level for each
> census tract. So a tesselation of my window, really.  Through much
> wrangling (possibly poorly-done) I was able to turn each predictor into
> a funxy--therefore they are essentially step functions, constant across
> a census tract and with abrupt changes at census tract boundaries.  I
> notice that rhohat calls for a continuous covariate. Could that be an
> issue?

I don't think so.  The rhohat() function works by smoothing, and
smoothing a step function doesn't really make sense.  But I think that
rhohat() should still give you an answer, even though it's a crappy answer.

> Although, I have one predictor that is a continuous distfun, and
> I get the same error message from rhohat with that.

Don't think that I can come up with any useful suggestions as to how to
proceed from here.  Sorry.

cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Fri, 12/15/2017 - 16:43
Thanks Rolf. I'm going to have to reflect more on my code and my data,
to understand better what is going on.

Obviously this won't help you much, without having access to all my data
and preceeding code, but the error message that is tripping me up is:

> rhohat(m12, pov.f)
Error:  the fitted intensity contains NA values

And yet,

table(is.na(fitted(m12)))
FALSE
  876

The predicted intensity, however, contains many NA values:

table(is.na(predict(m12)$v))
FALSE  TRUE
10379  6005

I try to force predictions only within my window by specifying locations
(which I think requires a binary mask window) but get the same result:

> table(is.na(predict(m12, locations = as.mask(sremsWindow))$v))
FALSE  TRUE
10379  6005

Does rhohat use fitted values (at quadrature points) or predicted values
(on a 128 x 128 pixel grid within the window)? Top of p. 415 in your
book Spatial Point Patterns: Methodology and Applications wtih R seems
to indicate the latter, while the error message from my rhohat command
above speaks of fitted values.  And how is a rectangular 128 x 128 grid
fit in an irregularly-shaped polygonal window?  Maybe that's how NA
predicted values arise--pixels outside an intra-window rectangular grid
but still inside the window?


And I can see now that no residuals from the model are missing:

> table(is.na(residuals(m12)$val))
FALSE
  876

All the NA's in my predicted values *around* my window, but inside the
bounding rectangle, led me down the garden path.

The origin of most of my predictors, such as pov.f, are shapefiles from
the US Census Bureau, with a discrete value of poverty level for each
census tract. So a tesselation of my window, really.  Through much
wrangling (possibly poorly-done) I was able to turn each predictor into
a funxy--therefore they are essentially step functions, constant across
a census tract and with abrupt changes at census tract boundaries.  I
notice that rhohat calls for a continuous covariate. Could that be an
issue?  Although, I have one predictor that is a continuous distfun, and
I get the same error message from rhohat with that.


Thanks.

--Chris Ryan

Rolf Turner wrote:
> set.seed(42)
> X <- runifpoint(20,win=test.window)
> xxx <- ppm(X ~ test.im)
> plot(residuals(xxx))
> # No sign of any missing values.

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Fri, 12/15/2017 - 15:20

I am probably misunderstanding something (as usual) but I cannot fathom
what you *expect* the values of the funxy object or the image to *be*,
outside of the window.

See inline below.

On 16/12/17 08:37, Christopher W. Ryan wrote:
> Using R 3.3.3 and spatstat 1.50-0 on Windows 7.  MWE below.

Thank you for providing a clear and simple example.

>
> I have a SpatialPolygonsDataFrame of census tract poverty levels in 3
> contiguous counties in the US, called sremsPoverty. I want to use this
> as a predictor in a ppm model. The window for the point pattern is the
> three counties--so an irregular polygon--called sremsWindow.
>
> I understand ppm predictors need to be an image, a tesselation, a funxy,
> a window, or a single number. So I proceed as follows:
>
> ### Poverty
> p <- slot(sremsPoverty, "polygons")
> v <- lapply(p, function(z) { SpatialPolygons(list(z)) })
> sat <- tess(tiles = lapply(v, as.owin) )
> pov.f <- as.function.tess(sat, values = sremsPoverty@data$propPoverty)
>
> Thus pov.f is a spatial function (funxy) that I can use in a call to ppm():
>
> m1 <- ppm(unique.cases.ppp ~ pov.f)
>
> pov.f looks as expected when I plot it. But examing the structure of
> as.im(pov.f) it appears it is a 128 x 128 pixel array, with the value of
> the function at all pixels outside of the irregular polygonal window,
> but inside the bounding rectangle, set to NA.
What else could/should they be set to?

> I think this is the cause
> of the NA values I am seeing among the residuals from m1,

I don't think this is the case.  Perhaps more detail is needed; perhaps
Adrian or Ege will be able to provide more insight.

> and those NA
> residuals are causing me some trouble with model diagnostics such as
> rhohat().

Again I, at least, would need more detail before being able to provide
any constructive comment.

> How do I constrain the funxy (or the image I can derive from it) to the
> irregular polygonal window, so as to eliminate the NA values outside the
> window but inside the bounding rectangle? Or can I constrain the
> modeling activity of ppm() to the window?

The "modelling activity of ppm()" is *ALWAYS* constrained to the window
of the pattern to which ppm() is applied.  This is fundamental to the
the way that ppm() works, and to what a window *means*.

>
> Thanks.
>
> --Chris Ryan
> Broome County Health Department
> Binghamton University
> SUNY Upstate Medical University
> Binghamt, NY, US
>
> MWE:
>
> x <- c(0, 2.6, 3, 1, 0)
> y <- c(1,2,3,2,1)
> test.window <- owin(poly=list(x=x,y=y))
> plot(test.window)  ## looks as expected
> ## make spatial function
> test.f <- function(x,y){x+y}
> test.funxy <- funxy(test.f, W = test.window)  ## I *thought* this would
> restrict the funxy to the window, but I don't think it does
> plot(test.funxy)  ## looks as expected
> ## but the image from test.funxy is square, with NA values outside the
> polygonal window, which is not square
> test.im <- as.im(test.funxy)
> str(test.im)
Again, what could the values of the image, outside of test.window,
possibly *be*, other than NA?

> ## my incorrect attempts to restrict the image to the window yields a
> numeric vector
> str(test.im[test.window])

You need to do

      new.im <- test.im[test.window,drop=FALSE]

to get an image rather than a numeric vector.  However the "new.im" that
you obtain will be identical to test.im:

     > all.equal(new.im,test.im)
     [1] TRUE

> ## or an error message
> window(test.im) <- test.window

You need a capital "W" on Window(test.im) (to avoid an error message).
But again this won't change anything.

Finally:

set.seed(42)
X <- runifpoint(20,win=test.window)
xxx <- ppm(X ~ test.im)
plot(residuals(xxx))
# No sign of any missing values.

cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

how to limit an image or a funxy to an irregular polygonal window, instead of the whole enclosing rectangle?

Fri, 12/15/2017 - 13:37
Using R 3.3.3 and spatstat 1.50-0 on Windows 7.  MWE below.

I have a SpatialPolygonsDataFrame of census tract poverty levels in 3
contiguous counties in the US, called sremsPoverty. I want to use this
as a predictor in a ppm model. The window for the point pattern is the
three counties--so an irregular polygon--called sremsWindow.

I understand ppm predictors need to be an image, a tesselation, a funxy,
a window, or a single number. So I proceed as follows:

### Poverty
p <- slot(sremsPoverty, "polygons")
v <- lapply(p, function(z) { SpatialPolygons(list(z)) })
sat <- tess(tiles = lapply(v, as.owin) )
pov.f <- as.function.tess(sat, values = sremsPoverty@data$propPoverty)

Thus pov.f is a spatial function (funxy) that I can use in a call to ppm():

m1 <- ppm(unique.cases.ppp ~ pov.f)

pov.f looks as expected when I plot it. But examing the structure of
as.im(pov.f) it appears it is a 128 x 128 pixel array, with the value of
the function at all pixels outside of the irregular polygonal window,
but inside the bounding rectangle, set to NA. I think this is the cause
of the NA values I am seeing among the residuals from m1, and those NA
residuals are causing me some trouble with model diagnostics such as
rhohat().

How do I constrain the funxy (or the image I can derive from it) to the
irregular polygonal window, so as to eliminate the NA values outside the
window but inside the bounding rectangle? Or can I constrain the
modeling activity of ppm() to the window?

Thanks.

--Chris Ryan
Broome County Health Department
Binghamton University
SUNY Upstate Medical University
Binghamt, NY, US

MWE:

x <- c(0, 2.6, 3, 1, 0)
y <- c(1,2,3,2,1)
test.window <- owin(poly=list(x=x,y=y))
plot(test.window)  ## looks as expected
## make spatial function
test.f <- function(x,y){x+y}
test.funxy <- funxy(test.f, W = test.window)  ## I *thought* this would
restrict the funxy to the window, but I don't think it does
plot(test.funxy)  ## looks as expected
## but the image from test.funxy is square, with NA values outside the
polygonal window, which is not square
test.im <- as.im(test.funxy)
str(test.im)

## my incorrect attempts to restrict the image to the window yields a
numeric vector
str(test.im[test.window])
## or an error message
window(test.im) <- test.window

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Intensity Estimation Methods

Fri, 12/15/2017 - 06:23
Thanks a lot.
I did superimpose the original pattern to pixel images of intensities. With adaptive smothing intensity higher zones are too narrowed.
Also the differences of intensities are getting higher in the study region. The fixed bandwidth choosen with bw.ppl( ) in spatstat give me a better Picture. This is my opinion and also this is in my case . I use this bandwidth as a global bandwidth for adaptive smoothing.


-----Original Message-----
From: Rolf Turner [mailto:[hidden email]]
Sent: Thursday, December 14, 2017 11:14 PM
To: Cenk İÇÖZ <[hidden email]>
Cc: [hidden email]; [hidden email]; Ege Rubak <[hidden email]>
Subject: Re: [R-sig-Geo] Intensity Estimation Methods

On 15/12/17 01:01, Cenk İÇÖZ via R-sig-Geo wrote:
> Hello list,,
>
> I have a spatial point pattern . I am trying to estimate its intensity
> both with a fixed bandwidth and  with an adaptive bandwidth. How could
> I compare the goodness of these two fits? I mean are there any things
> like mse, aic or any other criteria??? I want to compare the
> difference between the estimated intensity and the original pattern's
> intensity.

I think that the following fortune (fortunes::fortune(340)) might be
relevant:

> Bandwidth selection is an unresolved (and possibly unsolvable) problem
> in smoothing, so you're perfectly justified in trying/choosing an
> arbitrary value if it produces good pictures!
>    -- Adrian Baddeley (answering a user's question about the choice of
>       smoothing parameter when using the density.ppp() function from the
>       spatstat package)
>       private communication (March 2013)

I am cc-ing to the man himself to see if he has further comment.

cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276


________________________________

Bu elektronik posta ve onunla iletilen bütün dosyalar sadece yukarıda isimleri belirtilen kişiler arasında özel haberleşme amacını taşımakta olup gönderici tarafından alınması amaçlanan yetkili gerçek ya da tüzel kişinin kullanımına aittir. Eğer bu elektronik posta size yanlışlıkla ulaşmışsa, elektronik postanın içeriğini açıklamanız, kopyalamanız, yönlendirmeniz ve kullanmanız kesinlikle yasaktır. Bu durumda, lütfen mesajı geri gönderiniz ve sisteminizden siliniz. Anadolu Üniversitesi bu mesajın içerdiği bilgilerin doğruluğu veya eksiksiz olduğu konusunda herhangi bir garanti vermemektedir. Bu nedenle bu bilgilerin ne şekilde olursa olsun içeriğinden, iletilmesinden, alınmasından ve saklanmasından sorumlu değildir. Bu mesajdaki görüşler yalnızca gönderen kişiye aittir ve Anadolu Üniversitesinin görüşlerini yansıtmayabilir.

This electronic mail and any files transmitted with it are intended for the private use of the people named above. If you are not the intended recipient and received this message in error, forwarding, copying or use of any of the information is strictly prohibited. Any dissemination or use of this information by a person other than the intended recipient is unauthorized and may be illegal. In this case, please immediately notify the sender and delete it from your system. Anadolu University does not guarantee the accuracy or completeness of any information included in this message. Therefore, by any means Anadolu University is not responsible for the content of the message, and the transmission, reception, storage, and use of the information. The opinions expressed in this message only belong to the sender of it and may not reflect the opinions of Anadolu University.

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Intensity Estimation Methods

Thu, 12/14/2017 - 14:13
On 15/12/17 01:01, Cenk İÇÖZ via R-sig-Geo wrote:
> Hello list,,
>
> I have a spatial point pattern . I am trying to estimate its
> intensity both with a fixed bandwidth and  with an adaptive
> bandwidth. How could I compare the goodness of these two fits? I mean
> are there any things like mse, aic or any other criteria??? I want to
> compare the difference between the estimated intensity and the
> original pattern's intensity.

I think that the following fortune (fortunes::fortune(340)) might be
relevant:

> Bandwidth selection is an unresolved (and possibly unsolvable) problem in
> smoothing, so you're perfectly justified in trying/choosing an arbitrary value
> if it produces good pictures!
>    -- Adrian Baddeley (answering a user's question about the choice of
>       smoothing parameter when using the density.ppp() function from the
>       spatstat package)
>       private communication (March 2013)

I am cc-ing to the man himself to see if he has further comment.

cheers,

Rolf Turner

--
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: comparing one raster to a stack and condition

Thu, 12/14/2017 - 08:20
Thanks Barry,
 
this is ingenious ! It works really well. For the record below is my code
using this and writing the result to disk.
Karsten
 
 
library(raster)
# Use pattern arg to return a wildcard to get a list of all tifs in dir
alltiffs = list.files(getwd(), pattern="*\\.tif$", full.names=TRUE)
januarygrids =  alltiffs[grepl("*01.tif*", alltiffs)]

# Create raster stack of january rain grids
r <- stack(januarygrids, quick=TRUE)
 
# set current january layer to compare with
currentmonth <- "es_af.2017.01.tif"
currentmonthraster <- raster(currentmonth)
 
# create count for each cell where january rain is below the rain value in
the raster stack layers
resultbelow01 <- sum(r<currentmonthraster)

# write resulttif
writeRaster(resultbelow01, filename='below01.tif', overwrite=TRUE)

  _____  

From: [hidden email] [mailto:[hidden email]] On Behalf Of
Barry Rowlingson
Sent: Donnerstag, 14. Dezember 2017 09:29
To: karsten
Cc: r-sig-geo
Subject: Re: [R-sig-Geo] comparing one raster to a stack and condition


Here's a way - first let's make some sample data in a stack:

 maker = function(d){raster(matrix(runif(16),4,4))}
 rains = stack(lapply(1:10, maker))


so `rains` is a stack of 10 4x4 rasters with random numbers in. Now the
raster we want to test:


 r1 = maker()


Okay, all set up. We have a raster and a stack, then:

 r1 < rains

is a stack of 1s and 0s where r1 is less than the cell in each layer of
rains, and then we can do:

 sum(r1 < rains)

to total those up.

plot(sum(r1<rains))


should map how many times r1 is less than the values in rains.

Barry



On Wed, Dec 13, 2017 at 9:13 PM, karsten <[hidden email]> wrote:


Hi All,

I am trying to compare one precipitation raster to a stack of precipitation
raster and would like to create a result raster with a count of how often
the raster value is below those of the stack.
So far I have the following code:

----------------------------------------------

library(zoo)
library(raster)
# create raster stack for Januaries
alltiffs = list.files(getwd(), pattern="*\\.tif$", full.names=TRUE)

#filter the ones with 01 in it for january
januarygrids =  alltiffs[grepl("*.01.*", alltiffs)]

# Create raster stack of grids
r <- stack(januarygrids, quick=TRUE)

# set current january layer to compare with
currentmonth <- "es_af.2017.01.tif"
currentmonthraster <- raster(currentmonth)

# function to count how oftetn current ratser is below values in stack,
input r and currentmonthraster
belowCurrentRaster <- function(x, y) {
  sum(x > y)
}

-------------------------------------------------------

Now I thought I could use zApply or overlay to get the count from the
belowCurrentRaster function and write those counts into a result raster.
But I could not figure out how to make this work.
Any ideas appreciated.

Cheers
Karsten Vennemann
Terra GIS LTD


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/
<https://stat.ethz.ch/mailman/listinfo/r-sig-geo> listinfo/r-sig-geo




        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Intensity Estimation Methods

Thu, 12/14/2017 - 06:01
Hello list,,

I have a spatial point pattern . I am trying to estimate its intensity both with a fixed bandwidth and  with an adaptive bandwidth.
How could I compare the goodness of these two fits? I mean are there any things like mse, aic or any other criteria???
I want to compare the difference between the estimated intensity and the original pattern's intensity.

Thanks
Cenk





________________________________

Bu elektronik posta ve onunla iletilen bütün dosyalar sadece yukarıda isimleri belirtilen kişiler arasında özel haberleşme amacını taşımakta olup gönderici tarafından alınması amaçlanan yetkili gerçek ya da tüzel kişinin kullanımına aittir. Eğer bu elektronik posta size yanlışlıkla ulaşmışsa, elektronik postanın içeriğini açıklamanız, kopyalamanız, yönlendirmeniz ve kullanmanız kesinlikle yasaktır. Bu durumda, lütfen mesajı geri gönderiniz ve sisteminizden siliniz. Anadolu Üniversitesi bu mesajın içerdiği bilgilerin doğruluğu veya eksiksiz olduğu konusunda herhangi bir garanti vermemektedir. Bu nedenle bu bilgilerin ne şekilde olursa olsun içeriğinden, iletilmesinden, alınmasından ve saklanmasından sorumlu değildir. Bu mesajdaki görüşler yalnızca gönderen kişiye aittir ve Anadolu Üniversitesinin görüşlerini yansıtmayabilir.

This electronic mail and any files transmitted with it are intended for the private use of the people named above. If you are not the intended recipient and received this message in error, forwarding, copying or use of any of the information is strictly prohibited. Any dissemination or use of this information by a person other than the intended recipient is unauthorized and may be illegal. In this case, please immediately notify the sender and delete it from your system. Anadolu University does not guarantee the accuracy or completeness of any information included in this message. Therefore, by any means Anadolu University is not responsible for the content of the message, and the transmission, reception, storage, and use of the information. The opinions expressed in this message only belong to the sender of it and may not reflect the opinions of Anadolu University.

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: comparing one raster to a stack and condition

Thu, 12/14/2017 - 02:28
Here's a way - first let's make some sample data in a stack:

 maker = function(d){raster(matrix(runif(16),4,4))}
 rains = stack(lapply(1:10, maker))

so `rains` is a stack of 10 4x4 rasters with random numbers in. Now the
raster we want to test:

 r1 = maker()

Okay, all set up. We have a raster and a stack, then:

 r1 < rains

is a stack of 1s and 0s where r1 is less than the cell in each layer of
rains, and then we can do:

 sum(r1 < rains)

to total those up.

plot(sum(r1<rains))

should map how many times r1 is less than the values in rains.

Barry



On Wed, Dec 13, 2017 at 9:13 PM, karsten <[hidden email]> wrote:

> Hi All,
>
> I am trying to compare one precipitation raster to a stack of precipitation
> raster and would like to create a result raster with a count of how often
> the raster value is below those of the stack.
> So far I have the following code:
>
> ----------------------------------------------
>
> library(zoo)
> library(raster)
> # create raster stack for Januaries
> alltiffs = list.files(getwd(), pattern="*\\.tif$", full.names=TRUE)
>
> #filter the ones with 01 in it for january
> januarygrids =  alltiffs[grepl("*.01.*", alltiffs)]
>
> # Create raster stack of grids
> r <- stack(januarygrids, quick=TRUE)
>
> # set current january layer to compare with
> currentmonth <- "es_af.2017.01.tif"
> currentmonthraster <- raster(currentmonth)
>
> # function to count how oftetn current ratser is below values in stack,
> input r and currentmonthraster
> belowCurrentRaster <- function(x, y) {
>   sum(x > y)
> }
>
> -------------------------------------------------------
>
> Now I thought I could use zApply or overlay to get the count from the
> belowCurrentRaster function and write those counts into a result raster.
> But I could not figure out how to make this work.
> Any ideas appreciated.
>
> Cheers
> Karsten Vennemann
> Terra GIS LTD
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

comparing one raster to a stack and condition

Wed, 12/13/2017 - 15:13
Hi All,
 
I am trying to compare one precipitation raster to a stack of precipitation
raster and would like to create a result raster with a count of how often
the raster value is below those of the stack.
So far I have the following code:
 
----------------------------------------------
 
library(zoo)
library(raster)
# create raster stack for Januaries
alltiffs = list.files(getwd(), pattern="*\\.tif$", full.names=TRUE)

#filter the ones with 01 in it for january
januarygrids =  alltiffs[grepl("*.01.*", alltiffs)]
 
# Create raster stack of grids
r <- stack(januarygrids, quick=TRUE)
 
# set current january layer to compare with
currentmonth <- "es_af.2017.01.tif"
currentmonthraster <- raster(currentmonth)
 
# function to count how oftetn current ratser is below values in stack,
input r and currentmonthraster
belowCurrentRaster <- function(x, y) {
  sum(x > y)
}
 
-------------------------------------------------------
 
Now I thought I could use zApply or overlay to get the count from the
belowCurrentRaster function and write those counts into a result raster.
But I could not figure out how to make this work.
Any ideas appreciated.
 
Cheers
Karsten Vennemann
Terra GIS LTD


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: +towgs84 in st_write

Wed, 12/13/2017 - 08:09
Dear Roger and Edzer

Many thanks for looking into this. Projection is perfectly transferred from R to QGIS when using a geopackage, so this is definitively a shapefile issue. It's absolutely time to change ... and since st_write to GPKG works perfectly now (thanks to your efforts) there is no obstacle to do so.
All best wishes
Manuel

-----Ursprüngliche Nachricht-----
Von: R-sig-Geo [mailto:[hidden email]] Im Auftrag von Roger Bivand
Gesendet: Mittwoch, 13. Dezember 2017 07:09
An: Edzer Pebesma <[hidden email]>
Cc: [hidden email]
Betreff: Re: [R-sig-Geo] +towgs84 in st_write

On Tue, 12 Dec 2017, Edzer Pebesma wrote:

>
>
> On 12/12/2017 11:27 PM, Roger Bivand wrote:
>> Does sf use morphToESRI()?
>
> No.
>

I think the "ESRI Shapefile" driver changed - at one stage it was needed
during writing some time ago, I think. Now the driver simply does it
itself (line 763 in ogrsf_frmts/ogrshapedatasource.cpp:

* -------------------------------------------------------------------- */
/*      Create the .prj file, if required.                              */
/* -------------------------------------------------------------------- */
     if( poSRS != NULL )
     {
         CPLString osPrjFile =
             CPLFormFilename( NULL, pszFilenameWithoutExt, "prj");

         // The shape layer needs its own copy.
         poSRS = poSRS->Clone();
         poSRS->morphToESRI();

         char *pszWKT = NULL;
         VSILFILE *fp = NULL;
         if( poSRS->exportToWkt( &pszWKT ) == OGRERR_NONE
             && (fp = VSIFOpenL( osPrjFile, "wt" )) != NULL )
         {
             VSIFWriteL( pszWKT, strlen(pszWKT), 1, fp );
             VSIFCloseL( fp );
         }

         CPLFree( pszWKT );

         poSRS->morphFromESRI();
     }

So the driver is doing what it believes ArcGIS would like to read - the
*.prj file isn't well specified. In https://issues.qgis.org/issues/2154 
Frank Warmerdam wrote eight years ago: "Yes, OGR does strip the towgs84
parameter when writing the .prj file. TOWGS84 is not a legal construct in
an ESRI Projection Engine string (for .prj files)."

Sounds like another reason to abandon shapefiles as not fit for purpose.

Roger

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: +towgs84 in st_write

Wed, 12/13/2017 - 00:09
On Tue, 12 Dec 2017, Edzer Pebesma wrote:

>
>
> On 12/12/2017 11:27 PM, Roger Bivand wrote:
>> Does sf use morphToESRI()?
>
> No.
>

I think the "ESRI Shapefile" driver changed - at one stage it was needed
during writing some time ago, I think. Now the driver simply does it
itself (line 763 in ogrsf_frmts/ogrshapedatasource.cpp:

* -------------------------------------------------------------------- */
/*      Create the .prj file, if required.                              */
/* -------------------------------------------------------------------- */
     if( poSRS != NULL )
     {
         CPLString osPrjFile =
             CPLFormFilename( NULL, pszFilenameWithoutExt, "prj");

         // The shape layer needs its own copy.
         poSRS = poSRS->Clone();
         poSRS->morphToESRI();

         char *pszWKT = NULL;
         VSILFILE *fp = NULL;
         if( poSRS->exportToWkt( &pszWKT ) == OGRERR_NONE
             && (fp = VSIFOpenL( osPrjFile, "wt" )) != NULL )
         {
             VSIFWriteL( pszWKT, strlen(pszWKT), 1, fp );
             VSIFCloseL( fp );
         }

         CPLFree( pszWKT );

         poSRS->morphFromESRI();
     }

So the driver is doing what it believes ArcGIS would like to read - the
*.prj file isn't well specified. In https://issues.qgis.org/issues/2154 
Frank Warmerdam wrote eight years ago: "Yes, OGR does strip the towgs84
parameter when writing the .prj file. TOWGS84 is not a legal construct in
an ESRI Projection Engine string (for .prj files)."

Sounds like another reason to abandon shapefiles as not fit for purpose.

Roger

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: +towgs84 in st_write

Tue, 12/12/2017 - 16:55


On 12/12/2017 11:27 PM, Roger Bivand wrote:
> Does sf use morphToESRI()?

No.
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: +towgs84 in st_write

Tue, 12/12/2017 - 16:27
On Tue, 12 Dec 2017, Edzer Pebesma wrote:

>
>
> On 12/12/2017 03:16 PM, [hidden email] wrote:
>> Dear list
>>
>> I have a feature projected in EPSG:2056 (http://spatialreference.org/ref/epsg/2056/). When st_write this to a shapefile I get a prj-File without the +towgs84 parameters and this results in an offset if the shapefile is projected on the fly in a QGIS project with OpenLayers (EPSG:3857). Placement is correct if I open in a project with EPSG:2056 or if I manually assign EPSG:2056 to the shapefile.
>> This looks a bit similar to this older thread https://www.mail-archive.com/r-sig-geo@.../msg06462.html but in this case the +towgs84 parameters are well defined.
>> The feature has
>>> st_crs(p)
>> Coordinate Reference System:
>>   EPSG: 2056
>>   proj4string: "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"
>> st_write generates a prj reading PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",GEOGCS["GCS_Bessel 1841",DATUM["D_unknown",SPHEROID["bessel",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],UNIT["Meter",1]]
>> The result is that QGIS interprets this as a user defined CRS with +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs . The +towgs84 are assumed to be 0, I guess, and this results in the offset.
>>
>> Does anybody know if there is a way to specify the information written by st_write to the .prj-File?
>
> I can reproduce that, and see the same when writing to ESRI Shapefile
> using rgdal::writeOGR. The funny thing is that from gdal one also gets
> the following back, when asking for proj.4 -> wkt conversion:
I think this is the writeOGR() morphToESRI argument (but commented out,
maybe in the "ESRI Shapefile" driver?). Does sf use morphToESRI()? This
changes the WKT SRS representation:

library(rgdal)
CRSargs(CRS("+init=epsg:2056"))
showWKT(CRSargs(CRS("+init=epsg:2056")), morphToESRI=TRUE)
showWKT(CRSargs(CRS("+init=epsg:2056")), morphToESRI=FALSE)

Roger

>
> cat(st_as_text(st_crs(2056), pretty = TRUE))
> PROJCS["unnamed",
>    GEOGCS["Bessel 1841",
>        DATUM["unknown",
>            SPHEROID["bessel",6377397.155,299.1528128],
>            TOWGS84[674.374,15.056,405.346,0,0,0,0]],
>        PRIMEM["Greenwich",0],
>        UNIT["degree",0.0174532925199433]],
>    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
>    PARAMETER["latitude_of_center",46.95240555555556],
>    PARAMETER["longitude_of_center",7.439583333333333],
>    PARAMETER["azimuth",90],
>    PARAMETER["rectified_grid_angle",90],
>    PARAMETER["scale_factor",1],
>    PARAMETER["false_easting",2600000],
>    PARAMETER["false_northing",1200000],
>    UNIT["Meter",1]]
>
> which has the towgs parameters.
>
> If I write this file to GeoPackage, i.e. by
>
> st_write(nc, "xx.gpkg")
>
> I do see the towgs parameters when looking at it with ogrinfo (or
> reading it back with st_read).
>
> One more reason to ditch shapefiles?
>
>>
>> Many thanks for suggestions
>> Manuel
>>
>>
>>
>> -----
>> Manuel Schneider (Ing.-Agr. ETH, Dr. sc. nat.)
>> Team leader, Mountain grassland ecology & management
>>
>> Federal Department of Economic Affairs, Education and Research EAER
>> Agroscope
>>
>> Reckenholzstrasse 191, CH-8046 Z???rich
>> Ph. +41 58 468 75 98
>> Fax  +41 58 468 72 01
>> [hidden email]<mailto:[hidden email]>
>> https://www.agroscope.admin.ch/agroscope/en/home/topics/plant-production/forage-grassland-grazing-systems.html
>> http://scholar.google.com/citations?user=a0CE8xoAAAAJ&hl=de
>>
>> www.agroscope.ch<http://www.agroscope.ch/>  I good food, healthy environment
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: +towgs84 in st_write

Tue, 12/12/2017 - 16:03


On 12/12/2017 03:16 PM, [hidden email] wrote:
> Dear list
>
> I have a feature projected in EPSG:2056 (http://spatialreference.org/ref/epsg/2056/). When st_write this to a shapefile I get a prj-File without the +towgs84 parameters and this results in an offset if the shapefile is projected on the fly in a QGIS project with OpenLayers (EPSG:3857). Placement is correct if I open in a project with EPSG:2056 or if I manually assign EPSG:2056 to the shapefile.
> This looks a bit similar to this older thread https://www.mail-archive.com/r-sig-geo@.../msg06462.html but in this case the +towgs84 parameters are well defined.
> The feature has
>> st_crs(p)
> Coordinate Reference System:
>   EPSG: 2056
>   proj4string: "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"
> st_write generates a prj reading PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",GEOGCS["GCS_Bessel 1841",DATUM["D_unknown",SPHEROID["bessel",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],UNIT["Meter",1]]
> The result is that QGIS interprets this as a user defined CRS with +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs . The +towgs84 are assumed to be 0, I guess, and this results in the offset.
>
> Does anybody know if there is a way to specify the information written by st_write to the .prj-File?
I can reproduce that, and see the same when writing to ESRI Shapefile
using rgdal::writeOGR. The funny thing is that from gdal one also gets
the following back, when asking for proj.4 -> wkt conversion:

cat(st_as_text(st_crs(2056), pretty = TRUE))
PROJCS["unnamed",
    GEOGCS["Bessel 1841",
        DATUM["unknown",
            SPHEROID["bessel",6377397.155,299.1528128],
            TOWGS84[674.374,15.056,405.346,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
    PARAMETER["latitude_of_center",46.95240555555556],
    PARAMETER["longitude_of_center",7.439583333333333],
    PARAMETER["azimuth",90],
    PARAMETER["rectified_grid_angle",90],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",2600000],
    PARAMETER["false_northing",1200000],
    UNIT["Meter",1]]

which has the towgs parameters.

If I write this file to GeoPackage, i.e. by

st_write(nc, "xx.gpkg")

I do see the towgs parameters when looking at it with ogrinfo (or
reading it back with st_read).

One more reason to ditch shapefiles?

>
> Many thanks for suggestions
> Manuel
>
>
>
> -----
> Manuel Schneider (Ing.-Agr. ETH, Dr. sc. nat.)
> Team leader, Mountain grassland ecology & management
>
> Federal Department of Economic Affairs, Education and Research EAER
> Agroscope
>
> Reckenholzstrasse 191, CH-8046 Z�rich
> Ph. +41 58 468 75 98
> Fax  +41 58 468 72 01
> [hidden email]<mailto:[hidden email]>
> https://www.agroscope.admin.ch/agroscope/en/home/topics/plant-production/forage-grassland-grazing-systems.html
> http://scholar.google.com/citations?user=a0CE8xoAAAAJ&hl=de
>
> www.agroscope.ch<http://www.agroscope.ch/>  I good food, healthy environment
>
>
>
> [[alternative HTML version deleted]]
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

+towgs84 in st_write

Tue, 12/12/2017 - 08:16
Dear list

I have a feature projected in EPSG:2056 (http://spatialreference.org/ref/epsg/2056/). When st_write this to a shapefile I get a prj-File without the +towgs84 parameters and this results in an offset if the shapefile is projected on the fly in a QGIS project with OpenLayers (EPSG:3857). Placement is correct if I open in a project with EPSG:2056 or if I manually assign EPSG:2056 to the shapefile.
This looks a bit similar to this older thread https://www.mail-archive.com/r-sig-geo@.../msg06462.html but in this case the +towgs84 parameters are well defined.
The feature has
> st_crs(p)
Coordinate Reference System:
  EPSG: 2056
  proj4string: "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"
st_write generates a prj reading PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",GEOGCS["GCS_Bessel 1841",DATUM["D_unknown",SPHEROID["bessel",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER["azimuth",90],PARAMETER["scale_factor",1],PARAMETER["false_easting",2600000],PARAMETER["false_northing",1200000],UNIT["Meter",1]]
The result is that QGIS interprets this as a user defined CRS with +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs . The +towgs84 are assumed to be 0, I guess, and this results in the offset.

Does anybody know if there is a way to specify the information written by st_write to the .prj-File?

Many thanks for suggestions
Manuel



-----
Manuel Schneider (Ing.-Agr. ETH, Dr. sc. nat.)
Team leader, Mountain grassland ecology & management

Federal Department of Economic Affairs, Education and Research EAER
Agroscope

Reckenholzstrasse 191, CH-8046 Z�rich
Ph. +41 58 468 75 98
Fax  +41 58 468 72 01
[hidden email]<mailto:[hidden email]>
https://www.agroscope.admin.ch/agroscope/en/home/topics/plant-production/forage-grassland-grazing-systems.html
http://scholar.google.com/citations?user=a0CE8xoAAAAJ&hl=de

www.agroscope.ch<http://www.agroscope.ch/>  I good food, healthy environment



        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

3D point pattern

Tue, 12/12/2017 - 06:01
Dear members,

I write in order to ask if someone could indicate me if there is a way to measure and compare between groups the 3D homogeneity of a distribution of points in 3D. I think to some extension in 3D of usual point pattern analysis.

Thanks in advance for any advice

Best

Paolo

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Pages