[R-sig-Geo] spatial clusters

Roger Bivand Roger.Bivand at nhh.no
Tue Dec 21 09:40:35 CET 2010


On Mon, 20 Dec 2010, dorina.lazar wrote:

>
> Dear all,
>
> Finally, with your help, I make some progress to apply SKATER algorithm to
> European countries.
>
> Here is the R code, following the example from  spdep package pdf:
>
>> country<-wrld_simpl[wrld_simpl$ISO3%in%c("AZE", "ALB", "ARM", "BIH",
> "BGR", "CYP", "DNK", "IRL", "EST", "AUT", "CZE", "FIN",  "FRA", "GEO",
> "DEU", "GRC", "HRV", "HUN", "ISL", "ITA", "LVA", "BLR", "LTU", "SVK",
> "LIE", "MLT", "BEL", "AND",  "LUX", "MCO", "NLD", "NOR", "POL", "PRT",
> "ROU", "MDA", "RUS", "SVN"),]
>> ske<-spCbind(country, rclasif)
>> ske at data
>> de <- data.frame(scale(ske at data[13]))   #col 13 contains the health
> indicator (from rclasif)
>> nbe <- poly2nb(country)
>> ce <- nbcosts(nbe, de)
>> le <- nb2listw(nbe, ce, style="B", zero.policy=TRUE)
> Error in nb2listw(nbe, ce, style = "B", zero.policy = TRUE) :
> neighbours and glist do not conform

Using your corrected list of countries, I see two/three problems:

1) some proximate polygons do not touch because of line generalisation in 
the simplification of the country polygons, use:

nbe <- poly2nb(country, snap=0.1)

to force neighbours to touch.

2) you are left with three real islands - Iceland, Cyprus, and Malta. To 
join them in, consider:

knn4 <- knn2nb(knearneigh(crds, 4, longlat=TRUE),
   row.names=row.names(country))
knn4s <- make.sym.nb(knn4)
nbe1 <- union.nb(nbe, knn4s)

3) nbcosts() does not like single-column data frames, as it is supposed to 
be for multivariate data, so:

de <- data.frame(rnorm(46))
ce <- nbcosts(nbe1, de)
which(sapply(ce, function(x) isTRUE(any(is.na(x)))))

gives an NA at the end of each line of cost distances, but:

de <- data.frame(rnorm(46), rnorm(46))
ce <- nbcosts(nbe1, de)
which(sapply(ce, function(x) isTRUE(any(is.na(x)))))
le <- nb2listw(nbe1, ce, style="B")

doesn't, nor does:

de <- data.frame(rnorm(46), rep(1, 46))
ce <- nbcosts(nbe1, de)
which(sapply(ce, function(x) isTRUE(any(is.na(x)))))
le <- nb2listw(nbe1, ce, style="B")

From there,

mst.bh <- mstree(le,5)
res1 <- skater(mst.bh[,1:2], de, 2)
plot(country, col=heat.colors(3)[res1$groups])

seem to work.

I think that it makes sense for nbcosts() to fail for non-connected 
observations, but perhaps it should explain why.

Hope this helps,

Roger


>
> I mention that:
> -	at "nbe", appears the message:  6 regions with no links: CYP DNK HRV ISL
> MLT SWE
> -	"ce" contain a lot of NA
> If you have some ideas about the cause of this Error, please help me to
> finalize my exercise on SKATER algoritm.
>
> Thanks a lot for your help,
> Dorina
>
>
>
>
> On Sat, 18 Dec 2010 18:42:35 +0100 (CET), Roger Bivand
> <Roger.Bivand at nhh.no> wrote:
>> On Sat, 18 Dec 2010, Roman Luštrik wrote:
>>
>>> Perhaps you could follow the instructions from this
>>>
> exchange<http://r-sig-geo.2731867.n2.nabble.com/Guam-island-map-td4969467.html#a4972725>and
>>> see if you can extract the desired country (assuming the resolution
>>> suits you).
>>>
>>
>> I don't think that we know which countries Dorina needs, or at which
>> administrative level. If there are only 30 country observations, and the
>
>> countries are not contiguous, then the whole exercise seems unnecessary.
>
>> Can Dorina use the simplified world map in the maptools package:
>>
>> library(maptools)
>> ?wrld_simpl
>>
>> adding the data using the ISO country names? If these are countries at
>> some administrative division level, then maybe www.gadm.org will have
> what
>> is needed, but the countries will have to be joined together, or the
> whole
>> map subsetted. If these are say EU NUTS2 regions, then use can be made
> of
>> GISCO shapefiles from Eurostat. But we don't know what boundaries are
>> needed, and of course to use Skater for clustering, you need the
>> boundaries.
>>
>> Roger
>>
>>>
>>> Cheers,
>>> Roman
>>>
>>>
>>>
>>> On Sat, Dec 18, 2010 at 6:20 PM, Dorina Lazar <drnlazar at yahoo.com>
> wrote:
>>>
>>>> Georg,
>>>>
>>>> Thanks for being very patient with me.
>>>>
>>>> I have not a shapefile for countries.
>>>>
>>>> I have tried to obtain it. I have loaded a shapefile for Finland, GIS
>>>> data,
>>>> from
>>>>
> http://www.eea.europa.eu/data-and-maps/data/eea-reference-grids/zipped-shape-file-finland
>>>>
>>>> and I tried to read the file FI_1K.shp
>>>> Fin<- readShapePoly(system.file("etc/shapes/ FI_1K.shp",
>>>> package="spdep"))
>>>> but...
>>>> Error in getinfo.shape(filen) : Error opening SHP file
>>>>
>>>> How can I obtain, for each country,  the shapefile spatial polygon in
> R
>>>> (SpatialPolygons-DataFrame)?
>>>>
>>>> Thanks again,
>>>> Dorina
>>>>
>>>>
>>>> --- On Thu, 12/16/10, Georg Ru? <research at georgruss.de> wrote:
>>>>
>>>>> From: Georg Ru? <research at georgruss.de>
>>>>> Subject: Re: [R-sig-Geo] spatial clusters
>>>>> To: "Dorina Lazar" <dorina.lazar at econ.ubbcluj.ro>
>>>>> Cc: "R-sig-geo mailing list" <r-sig-geo at r-project.org>
>>>>> Date: Thursday, December 16, 2010, 10:41 AM
>>>>> On 16/12/10 09:06:33, dorina.lazar
>>>>> wrote:
>>>>>> How to create a shapefile including my data about some
>>>>> socio-demographic
>>>>>> indicators (7),  for about 30 countries around
>>>>> the world (to be read with
>>>>>> readShapePoly)?
>>>>>
>>>>> Hi Dorina,
>>>>>
>>>>> I'm not quite sure what you have and what you want.
>>>>>
>>>>> I assume that you'd like to create a shapefile such that
>>>>> for every country
>>>>> (which is represented as a polygon) there's a vector of
>>>>> seven attributes
>>>>> attached to it:
>>>>>
>>>>> country1:  attr1 attr2 attr3 attr4 attr5 attr6 attr7
>>>>> country2:  attr1 attr2 attr3 attr4 attr5 attr6 attr7
>>>>> ...
>>>>>
>>>>> Do you have your countries as spatial polygons in R
>>>>> (SpatialPolygons-
>>>>> DataFrame or similar)? If yes, you should be able to add
>>>>> your socio-demo-
>>>>> graphic indicators to this data frame (just as in normal
>>>>> data frames). In
>>>>> this way you'd join the data you have inside R.
>>>>>
>>>>> If you have the shapefile for the countries available and
>>>>> would like to
>>>>> edit this directly to add the indicators, I'd suggest using
>>>>> a GIS like
>>>>> GRASS http://grass.fbk.eu/ Afterwards, this can be easily
>>>>> read into R.
>>>>>
>>>>>> Would be useful to suggest some good introductory
>>>>> lectures in spatial
>>>>>> statistics (links or books).
>>>>>
>>>>> I think the following is what Roger would suggest:
>>>>>
>>>>>  author = {Bivand, Roger S. and Pebesma, Edzer J. and
>>>>> G?mez-Rubio, Virgilio},
>>>>>  title = {Applied Spatial Data Analysis with R},
>>>>>  series = {Use R},
>>>>>  publisher = {Springer},
>>>>>
>>>>> http://www.asdar-book.org/
>>>>> http://www.springerlink.com/content/uw07v1/
>>>>> (Depending on your login/institute you may have fulltext
>>>>> access, but the
>>>>> book is definitely worth buying anyway.)
>>>>>
>>>>> Regards,
>>>>> Georg.
>>>>> --
>>>>> Research Assistant
>>>>> Otto-von-Guericke-Universit?t Magdeburg
>>>>> research at georgruss.de
>>>>> http://research.georgruss.de
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at r-project.org
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>
>>>
>>>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no


More information about the R-sig-Geo mailing list