Re: [R-sig-Geo] removing internal boundaries from a spatstat window that was created from three adjacent counties

[Rolf Turner]

On 13/02/16 11:59, Christopher W. Ryan wrote:


In spatstat, I want to create a single window from three adjacent
counties in New York State, USA.  My original data is a shapefile,
"cty036" showing all the New York State county boundaries.  Here's what
I've done so far:

== begin code =
library(shapefiles)
library(spatstat)
library(maptools)
library(sp)
library(RColorBrewer)
library(rgdal)

nysArea <- readOGR("C:/DATA/GeographicData", "cty036")

nysAreaUTM <- spTransform(nysArea, CRS("+proj=utm +zone=18 +ellps=GRS80
+units=m +no_defs"))

sremsAreaUTM <- subset(nysAreaUTM, NAME=="Broome" | NAME=="Tioga" |
NAME=="Chenango")

sremsWindow <- as(sremsAreaUTM, "owin")
 end code =


This works pretty well, producing a 3-county owin object. But the
dividing lines between the counties are shown in the window, whenever I
plot it or plot subsequent point patterns that use the window.  In
essence, in plots it looks like 3 polygons instead of one big one. I'd
prefer not to have the inter-county boundaries be visible--I'd rather
have just one big polygon for the whole area. How can I remove them?  Or
should I create the window differently from the get-go?

Thanks.


I think it probable that your owin object has come out as a multiple 
polygon.  Look at length(sremsWindow$bdry) --- this is probably 3 (or more).


I would guess that the internal boundaries actually consist of *pairs*
of closely juxtaposed lines --- which look like single lines when plotted.

Have you read the spatstat vignette "shapefiles"?  That might give you 
some guidance as to how to proceed.


I would have thought the automatic repair process that spatstat now uses 
would fix up this problem.  What version of spatstat are you using?


You *might* be able to fix things up by doing (something like):

W1 <- owin(poly=sremsWindow$bdry[[1]])
W2 <- owin(poly=sremsWindow$bdry[[2]])
W3 <- owin(poly=sremsWindow$bdry[[3]])

W  <- union.owin(W1,W2,W3)

But if my guess about the internal boundaries actually being pairs of 
line is correct, this may not work.


It's hard to say without having your sremsWindow object to experiment on.

Or perhaps you need to extract the three counties separately as owin
objects and the apply union.owin() to the three results.  (Again, if my 
guess is correct, this may not work.)


Adrian or Ege may be able to propose a more efficient/effective solution.

cheers,

Rolf Turner

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

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] different projection transformation R and gdal commandline

[Chris Reudenbach]

Hi Dominik,

If you use the gdalUtils package  there is no significant difference in 
the results using CLI or  R:


library(gdalUtils)
gdaltransform(s_srs="+proj=longlat +datum=WGS84", t_srs="+proj=lcc 
+lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0 +y_0=0  
+ellps=sphere  +a=637 +b=637 +units=m 
+no_defs",coords=matrix(c(-112.25,33.0), ncol = 2))

 [,1]  [,2] [,3]
[1,] -1306676 -629522.50


If I use  spTransform I can reproduce your results:


loc <- c(-112.25,33.0)
loc <- data.frame(matrix(unlist(loc), nrow=1, 
byrow=T),stringsAsFactors=FALSE)

colnames(loc)<-c("lon","lat")
coordinates(loc)<- ~lon+lat
proj4string(loc)<- CRS("+proj=longlat +datum=WGS84 +no_defs")
loc<-spTransform(loc,CRS('+proj=lcc +lat_1=28 +lat_2=50 
+lat_0=39.70001220694445 +lon_0=-98 +x_0=0 +y_0=0 +ellps=sphere 
+a=637 +b=637 +units=m +no_defs'))

loc@coords
   lon lat
1 -1306471 -629424

I suggest to focus on the sptransform() function

cheers Chris

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.3 LTS

locale:
 [1] LC_CTYPE=de_DE.UTF-8   LC_NUMERIC=C LC_TIME=de_DE.UTF-8
LC_COLLATE=de_DE.UTF-8
 [5] LC_MONETARY=de_DE.UTF-8LC_MESSAGES=de_DE.UTF-8 
LC_PAPER=de_DE.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C 
LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C


attached base packages:
[1] stats graphics  grDevices utils datasets  methods base

other attached packages:
[1] gdalUtils_2.0.1.7 raster_2.5-2  sp_1.2-2

loaded via a namespace (and not attached):
 [1] rgdal_1.1-3   tools_3.2.3   Rcpp_0.12.3 R.methodsS3_1.7.0 
codetools_0.2-14  grid_3.2.3
 [7] iterators_1.0.8   foreach_1.4.3 R.utils_2.2.0 
R.oo_1.19.0   lattice_0.20-33





Am 15.02.2016 um 21:37 schrieb Dominik Schneider:

Hi,
I'm struggling to use a custom projection. I am seeing differences with
someone using python proj4 bindings and when I compared my R results with
my commandline results I got even more confused. the coordinate
transformation is different for the two different methods.

could someone explain to me which one is wrong and why?
Thanks

R:
e=c(-112.25,-104.125,33,43.75)
box=as(extent(e),'SpatialPolygons')
proj4string(box)='+proj=longlat +datum=WGS84'
pstring='+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98
+x_0=0 +y_0=0 +ellps=sphere +a=637 +b=637 +units=m +no_defs'
xmin(spTransform(box,CRS(pstring)))
## [1] -1306471
ymin(spTransform(box,CRS(pstring)))
## [1] -713442.3

commandline:
iMac:~ $ gdaltransform -s_srs "+proj=longlat +datum=WGS84" -t_srs
"+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0
+y_0
=0  +ellps=sphere  +a=637 +b=637 +units=m +no_defs"
-112.25 33.
-1306675.75472246 -629522.472322824 0

[[alternative HTML version deleted]]

___
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] different projection transformation R and gdal commandline

[Dominik Schneider]

Sorry, forgot to add the details. rgdal and gdal 1.11.3 were installed from
kyngchaos.com

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-apple-darwin14.5.0 (64-bit)
Running under: OS X 10.11.3 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats graphics  grDevices utils datasets  methods
base

other attached packages:
 [1] rasterVis_0.37  latticeExtra_0.6-26 lattice_0.20-33
broom_0.4.0
 [5] tidyr_0.3.1 dplyr_0.4.3 doMC_1.3.4
 iterators_1.0.8
 [9] foreach_1.4.3   ncdf4_1.15  gridExtra_2.0.0
spdep_0.5-92
[13] Matrix_1.2-3ipred_0.9-5 MASS_7.3-45
RColorBrewer_1.1-2
[17] rgdal_1.0-7 stringr_1.0.0   ggplot2_2.0.0   plyr_1.8.3

[21] reshape2_1.4.1  raster_2.5-2sp_1.2-1
 gstat_1.1-0
[25] ProjectTemplate_0.6

loaded via a namespace (and not attached):
 [1] zoo_1.7-12   splines_3.2.3colorspace_1.2-6 spacetime_1.1-5
 [5] survival_2.38-3  prodlim_1.5.7hexbin_1.27.1DBI_0.3.1
 [9] lava_1.4.1   munsell_0.4.2gtable_0.1.2 codetools_0.2-14
[13] coda_0.18-1  psych_1.5.8  labeling_0.3 class_7.3-14
[17] xts_0.9-7Rcpp_0.12.2  scales_0.3.0 FNN_1.1
[21] deldir_0.1-9 mnormt_1.5-3 digest_0.6.8 stringi_1.0-1
[25] grid_3.2.3   tools_3.2.3  LearnBayes_2.15  magrittr_1.5
[29] lazyeval_0.1.10  assertthat_0.1   R6_2.1.1 rpart_4.1-10
[33] boot_1.3-17  intervals_0.15.1 nnet_7.3-11  nlme_3.1-122


Dominik Schneider
c 518.956.3978


On Mon, Feb 15, 2016 at 1:37 PM, Dominik Schneider <
dominik.schnei...@colorado.edu> wrote:

> Hi,
> I'm struggling to use a custom projection. I am seeing differences with
> someone using python proj4 bindings and when I compared my R results with
> my commandline results I got even more confused. the coordinate
> transformation is different for the two different methods.
>
> could someone explain to me which one is wrong and why?
> Thanks
>
> R:
> e=c(-112.25,-104.125,33,43.75)
> box=as(extent(e),'SpatialPolygons')
> proj4string(box)='+proj=longlat +datum=WGS84'
> pstring='+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98
> +x_0=0 +y_0=0 +ellps=sphere +a=637 +b=637 +units=m +no_defs'
> xmin(spTransform(box,CRS(pstring)))
> ## [1] -1306471
> ymin(spTransform(box,CRS(pstring)))
> ## [1] -713442.3
>
> commandline:
> iMac:~ $ gdaltransform -s_srs "+proj=longlat +datum=WGS84" -t_srs
> "+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0
> +y_0
> =0  +ellps=sphere  +a=637 +b=637 +units=m +no_defs"
> -112.25 33.
> -1306675.75472246 -629522.472322824 0
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] different projection transformation R and gdal commandline

[Dominik Schneider]

Hi,
I'm struggling to use a custom projection. I am seeing differences with
someone using python proj4 bindings and when I compared my R results with
my commandline results I got even more confused. the coordinate
transformation is different for the two different methods.

could someone explain to me which one is wrong and why?
Thanks

R:
e=c(-112.25,-104.125,33,43.75)
box=as(extent(e),'SpatialPolygons')
proj4string(box)='+proj=longlat +datum=WGS84'
pstring='+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98
+x_0=0 +y_0=0 +ellps=sphere +a=637 +b=637 +units=m +no_defs'
xmin(spTransform(box,CRS(pstring)))
## [1] -1306471
ymin(spTransform(box,CRS(pstring)))
## [1] -713442.3

commandline:
iMac:~ $ gdaltransform -s_srs "+proj=longlat +datum=WGS84" -t_srs
"+proj=lcc +lat_1=28 +lat_2=50 +lat_0=39.70001220694445 +lon_0=-98 +x_0=0
+y_0
=0  +ellps=sphere  +a=637 +b=637 +units=m +no_defs"
-112.25 33.
-1306675.75472246 -629522.472322824 0

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Geographic/Spatial Clustering - pyCluster for R?

[Timothée Giraud]

Hello,
You could have a look at the ClustGeo package :
>ClustGeo: Clustering of Observations with Geographical Constraints
>Functions which allow to integrate geographical constraints in Ward 
hierarchical clustering. Geographical maps of typologies obtained can be 
displayed with the use of shapefiles.


___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] R-sig-Geo Digest, Vol 150, Issue 14

[Corey Sparks]

Hi Maryia,
What function are you using to fit the models? If you're using the fitting 
routines in spdep, (lagsarlm, errorsarlm, sacsarlm, spautolm) then the AIC 
prints for each model using the summary() function. You can also get a pseudo 
R^2 using the summary(fit, Nagelkerke=T) option

for example:
require(spdep)
data(oldcol)
COL.lag.eig <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,
nb2listw(COL.nb, style="W"), method="eigen", quiet=FALSE)
summary(COL.lag.eig, Nagelkerke=T)

shows both the AIC and the pseudo R^2

Hope this helps

Corey Sparks, PhD
Associate Professor
Graduate Advisor of Record
Department of Demography
The University of Texas at San Antonio
501 West Cesar E Chavez Blvd
corey.sparks 'at' utsa.edu
coreysparks.weebly.com

On Feb 15, 2016, at 5:00 AM, r-sig-geo-requ...@r-project.org wrote:

Date: Sun, 14 Feb 2016 20:58:42 -0600
From: Maryia Bakhtsiyarava 
To: R-sig-geo Mailing List 
Subject: [R-sig-Geo] AIC/R^2 in splm
Message-ID:

Content-Type: text/plain; charset="UTF-8"

Dear list members,

I am estimating a spatial lag model with time-period fixed effects using
package splm. I would like to obtain some goodness-of-fit measures for my
models but I cannot figure out how to do it. The traditional AIC extraction
function doesn't work for a an object of class "splm".

The only thing I can extract is the log likelihood, using which in theory I
can calculate AIC, but even in that case I am not sure about the degrees of
freedom to use in the calculation (do I count time dummies, lag and
intercept as parameters?). I tried df.residual(model) but I got NULL.

Is there another way to obtain AIC and/or R^2? I am sure people encountered
this problem before, so if you have any advice on how to obtain model
statistics, I would greatly appreciate it.

Thank you,
Maryia
--
Maryia Bakhtsiyarava
Graduate student
Department of Geography, Environment and Society
University of Minnesota, Twin Cities

Research Assistant
TerraPop Project
Minnesota Population Center

414 Social Sciences, 267 19th Ave S, Minneapolis, MN 55455

[[alternative HTML version deleted]]



--

Message: 3
Date: Mon, 15 Feb 2016 09:33:44 +0100
From: Tobias R?ttenauer 
To: 
Cc: R-sig-Geo@r-project.org
Subject: Re: [R-sig-Geo] LM test for spatial dependence with panel
data
Message-ID: <001501d167cb$950379b0$bf0a6d10$@sowi.uni-kl.de>
Content-Type: text/plain; charset="iso-8859-1"

Dear Roger,

thanks for your answer!


Dear list members,

I'm currently estimating a fixed effects panel model and I want to
control for spatial dependence. Thus, I also estimated two spatial
fe-models, one with a spatial error term and one with a spatial error
term and spatial lag variable. Both lambda and rho are highly
significant and the independent variables of the models differ
considerably.

However, I would like to control the model specifications with a
Lagrange Multiplier test. More precisely, I would like to test for
spatial lag dependence (while allowing possible spatial error
dependence) and for spatial error dependence (while allowing possible
spatial lag dependence) like described by Anselin/ Bera/ Florax/ Yoon
(1996), which is possible for cross-sectional data using lm.LMtests
(spdep). Is there a possibility to estimate a similar test for panel
data? If I'm not mistaking, the splm-implemented tests do not account
for
spatial lag dependence?


You may be able to test the pooled model only using lm() on the data and
taking a Kronecker product of your spatial weights by a time-dimensioned
identity matrix, and using the function for cross-sectional models, but I
guess
that this is not what you want.

I mean, one could argue that a cross-sectional spatial dependence also leads
to spatially correlated variations over time, but this is not necessarily
the case. So you are right, it would be the best to implement an appropriate
test (if feasible).

You'd need to identify an appropriate test in
the literature, and see whether it can be implemented (and whether then it
actually works). In many cases, significant lag and error processes
together
are related to missing variables, and to spatial processes not matching
the
spatial units - aggregates - you are using.

You mean the error-correlation is a result of one or more omitted variables
that operate independently from the spatial relationship I modeled in the
spatial weights matrix? Then, it may be an idea to compare different
distance/neighborhood measures? Or do you mean it is a problem of the
aggregates itself (like too large spatial units)?

Thanks again,
Tobi


Hope this helps,

Roger

Thank you very much in advance for your help!

Best,
Tobi

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


--
Roger Bivand
Department of 

[R-sig-Geo] slope of regression line in moran scatter plot

[Takahiro Yoshida]

Dear list members,

In moran scatter plot, I checked
[1] Slope of lm((spatial lagged variable) ~ (variable)): moran coefficient
(spdep:::moran.plot),
[2] Slope of lm((spatial lagged standardized variable) ~ (standardized
variable)): moran coefficient, and
[3] Slope of lm((standardized spatial lagged variable) ~ (standardized
variable)): correlation coefficient.
However, I think Anselin (1999) remark case [3] is moran coefficient. Is it
true?

Anselin, L. (1999) Interactive techniques and exploratory spatial data
analysis. In: Longley P.A., Goodchild M.F., Maguire D.J., Rhind D.W. (eds)
Geographic information system: Principles, techniques, management and
applications, 253–266, Wiley, New York.

Thanks and regards,

Takahiro

--
Takahiro Yoshida
Ph.D. student
Real Estate & Spatial Statistics Laboratory
University of Tsukuba, Japan
E-mail: yoshida.takah...@sk.tsukuba.ac.jp

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] removing internal boundaries from a spatstat window that was created from three adjacent counties

[Christopher W. Ryan]

In spatstat, I want to create a single window from three adjacent
counties in New York State, USA.  My original data is a shapefile,
"cty036" showing all the New York State county boundaries.  Here's what
I've done so far:

== begin code =
library(shapefiles)
library(spatstat)
library(maptools)
library(sp)
library(RColorBrewer)
library(rgdal)

nysArea <- readOGR("C:/DATA/GeographicData", "cty036")

nysAreaUTM <- spTransform(nysArea, CRS("+proj=utm +zone=18 +ellps=GRS80
+units=m +no_defs"))

sremsAreaUTM <- subset(nysAreaUTM, NAME=="Broome" | NAME=="Tioga" |
NAME=="Chenango")

sremsWindow <- as(sremsAreaUTM, "owin")
 end code =


This works pretty well, producing a 3-county owin object. But the
dividing lines between the counties are shown in the window, whenever I
plot it or plot subsequent point patterns that use the window.  In
essence, in plots it looks like 3 polygons instead of one big one. I'd
prefer not to have the inter-county boundaries be visible--I'd rather
have just one big polygon for the whole area. How can I remove them?  Or
should I create the window differently from the get-go?

Thanks.

--Chris
-- 
Christopher W. Ryan, MD, MS
cryanatbinghamtondotedu
https://www.linkedin.com/in/ryancw

Early success is a terrible teacher. You’re essentially being rewarded
for a lack of preparation, so when you find yourself in a situation
where you must prepare, you can’t do it. You don’t know how.
--Chris Hadfield, An Astronaut's Guide to Life on Earth

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] LM test for spatial dependence with panel data

[Tobias Rüttenauer]

Dear Roger,

thanks for your answer!

 
> > Dear list members,
> >
> > I'm currently estimating a fixed effects panel model and I want to
> > control for spatial dependence. Thus, I also estimated two spatial
> > fe-models, one with a spatial error term and one with a spatial error
> > term and spatial lag variable. Both lambda and rho are highly
> > significant and the independent variables of the models differ
> considerably.
> >
> > However, I would like to control the model specifications with a
> > Lagrange Multiplier test. More precisely, I would like to test for
> > spatial lag dependence (while allowing possible spatial error
> > dependence) and for spatial error dependence (while allowing possible
> > spatial lag dependence) like described by Anselin/ Bera/ Florax/ Yoon
> > (1996), which is possible for cross-sectional data using lm.LMtests
> > (spdep). Is there a possibility to estimate a similar test for panel
> > data? If I'm not mistaking, the splm-implemented tests do not account
for
> spatial lag dependence?
> >
> 
> You may be able to test the pooled model only using lm() on the data and
> taking a Kronecker product of your spatial weights by a time-dimensioned
> identity matrix, and using the function for cross-sectional models, but I
guess
> that this is not what you want. 

I mean, one could argue that a cross-sectional spatial dependence also leads
to spatially correlated variations over time, but this is not necessarily
the case. So you are right, it would be the best to implement an appropriate
test (if feasible).

> You'd need to identify an appropriate test in
> the literature, and see whether it can be implemented (and whether then it
> actually works). In many cases, significant lag and error processes
together
> are related to missing variables, and to spatial processes not matching
the
> spatial units - aggregates - you are using.

You mean the error-correlation is a result of one or more omitted variables
that operate independently from the spatial relationship I modeled in the
spatial weights matrix? Then, it may be an idea to compare different
distance/neighborhood measures? Or do you mean it is a problem of the
aggregates itself (like too large spatial units)? 

Thanks again,
Tobi

> 
> Hope this helps,
> 
> Roger
> 
> > Thank you very much in advance for your help!
> >
> > Best,
> > Tobi
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > 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; fax +47 55 95 91 00
> e-mail: roger.biv...@nhh.no
> http://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
> http://depsy.org/person/434412

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo