[R-sig-Geo] Problems With Cartography 2 - Help

[Andrés Peralta]

Thank you Chris,

I found a way to connect the three islands in a previous R-sig-Geo answer
from Roger Bivand (
https://stat.ethz.ch/pipermail/r-sig-geo/2013-July/018868.html). I followed
the code given in this answer and it seems to work fine (I plotted it and
it looks ok).

x.nb <- poly2nb(carto2012u2)
summary(x.nb)

### 1) fIND AREAS WITH NO NEIGHBORS
no_neighs <- which(card(x.nb) == 0) # returns set cardinality.

### To check that you get the same result using counts of graph components:

nc <- n.comp.nb(x.nb)
table(nc$comp.id)
no_neighs_comps <- which(unname(table(nc$comp.id)==1))
match(no_neighs_comps, nc$comp.id)

### 2) For each island find nearest neighbour

k1nb <- knn2nb(knearneigh(coordinates(carto2012u2), k=1,
                          longlat=!is.projected(carto2012u2)))

### 3) Assign a neighbour relationship between the island and its nearest
neighbour

is.symmetric.nb(x.nb, force=TRUE)
nb1 <- x.nb

res <- k1nb[no_neighs]
nb1[no_neighs] <- res
attr(nb1, "sym") <- is.symmetric.nb(nb1, force=TRUE)

# re-assign the now incorrect symmetry attribute

nb1

nb2 <- make.sym.nb(nb1)
is.symmetric.nb(nb2, force=TRUE)
nb2

table(n.comp.nb(nb2)$comp.id)

###

veins_inla <- nb2INLA("carto2012_nb.inla", nb2)

plot (carto2012u2, axes=F)
coords <- coordinates(carto2012u2)
plot(nb2, coords, add=T, col="red")

Hope this would be helpful for others with the same issue.


Andrés


On Thu, May 11, 2017 at 3:25 PM, chris english <
englishchristophera at gmail.com> wrote:

> Andres,
>
> I'm thinking you want to 'unsimplify' the topology prior to your poly2nb
> by a slight negative gBuffer on the Galapagos polygons,
> reduce the Galapagos polygons by the buffer and un-share the boundaries.
> I'm trying some code, but looking at a mix of these as
> guidance:
>
> https://www.rdocumentation.org/packages/rgeos/versions/0.
> 3-22/topics/gBuffer
> https://gist.github.com/mstrimas/1b4a4b93a9d4a158bce4
> https://gis.stackexchange.com/questions/93096/how-to-
> perform-a-true-gis-clip-of-polygons-layer-using-a-polygon-layer-in-r
>
> Saludos,
> Chris
>
> On Wed, May 10, 2017 at 4:27 AM, Andrés Peralta <tirico85 at gmail.com>
> wrote:
>
>> Hi to everyone,
>>
>> I`m working with the second administrative level cartography from Ecuador
>> (available in: http://www.ecuadorencifras.gob.ec//documentos/web-inec/
>> Cartografia/2015/Clasificador_Geografico/2012/SHP.zip
>> <http://www.ecuadorencifras.gob.ec//documentos/web-inec/Cartografia/2015/Clasificador_Geografico/2012/SHP.zip>).
>> The shape file is
>> called nxcantones.shp. I`ve been having a lot of problems opening and
>> working with this cartography in R; but i can open it in Q-GIS and ARCGIS
>> without any problem (I have opened it in both programs and saved it
>> again).
>> As the cartography has various objects with the same ID - aparently the
>> same areas but repeated; we had to run the following sintax in order to
>> have one ID in each area:
>>
>> *carto2012 <- readOGR("CANTONES2012/2012CLEAN.shp", "2012CLEAN",
>> stringsAsFactors=F) *
>> *proj4string(carto2012) <- CRS("+proj=utm +zone=17 +south +datum=WGS84
>> +units=m +no_defs")*
>>
>> *carto2012x <- unionSpatialPolygons(carto2012, IDs=carto2012$DPA_CANTON)*
>> *proj4string(carto2012x) <- CRS("+proj=utm +zone=17 +south +datum=WGS84
>> +units=m +no_defs")*
>>
>>
>> *datos <- data.frame(ID=row.names(carto2012x), stringsAsFactors = F)*
>> *row.names(datos) <- row.names(carto2012x)*
>> *carto2012x2 <- SpatialPolygonsDataFrame(carto2012x, data=datos)*
>>
>> *carto2012 <- carto2012x2 *
>>
>> For our work, we had to make a neighborhood matrix of the cartography.
>>
>> *x.nb <- poly2nb(carto2012u2) # U2 is the cartography after the recode of
>> some areas that had to be joined.*
>>
>> *summary(x.nb)*
>>
>> This seems to work fine. The problem is we need to connect the three
>> island
>> areas (the Galapagos) but when we try to do so, it connects many areas
>> along all the cartography. I've tried to do it manually using *edit.nb*
>> but
>> it does not work. I´ve also tried the following sintaxes:
>>
>> x.nb1 <- x.nb
>> which(card(x.nb1)==0) #to discover the id of these areas without
>> connections
>> id <- function(x){which(carto2012u2$ID==x)}
>>
>> x.nb1[[id(2001)]] = unique(as.integer(sort(c(x.nb1[[id(2001)]],
>> id(2002)))))
>> x.nb1[[id(2002)]] = unique(as.integer(sort(c(x.nb1[[id(2002)]],
>> id(2001)))))
>>
>> x.nb1[[id(2001)]] = unique(as.integer(sort(c(x.nb1[[id(2001)]],
>> id(2003)))))
>> x.nb1[[id(2003)]] = unique(as.integer(sort(c(x.nb1[[id(2003)]],
>> id(2001)))))
>>
>> x.nb1[[id(2002)]] = unique(as.integer(sort(c(x.nb1[[id(2002)]],
>> id(2003)))))
>> x.nb1[[id(2003)]] = unique(as.integer(sort(c(x.nb1[[id(2003)]],
>> id(2002)))))
>>
>> ####
>> x.nb1[[id("2001")]] = unique(as.integer(sort(c(x.nb1[[id("2001")]],
>> id("2002")))))
>> x.nb1[[id("2002")]] = unique(as.integer(sort( id("2001"))))
>>
>> x.nb1[[id("2001")]] = unique(as.integer(sort(c(x.nb1[[id("2001")]],
>> id("2003")))))
>> x.nb1[[id("2003")]] = unique(as.integer(sort( id("2001"))))
>>
>> x.nb1[[id("2002")]] = unique(as.integer(sort(c(x.nb1[[id("2002")]],
>> id("2003")))))
>> x.nb1[[id("2003")]] = unique(as.integer(sort( id("2002"))))
>>
>> Any ideas or suggestions? Your help will really be appreciated.
>>
>>
>> --
>>
>> *             Andrés Peralta*
>>         Pre-doctoral Researcher
>>      GREDS/EMCONET / ASPB
>> <https://www.upf.edu/greds-emconet/en/>   <http://aspb.cat/>
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>


-- 

*             Andrés Peralta*
        Pre-doctoral Researcher
     GREDS/EMCONET / ASPB
<https://www.upf.edu/greds-emconet/en/>   <http://aspb.cat/>

	[[alternative HTML version deleted]]