Problems With Cartography 2 - Help

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Problems With Cartography 2 - Help

Andrés Peralta
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). 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Problems With Cartography 2 - Help

chris english-2
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 <[hidden email]> 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). 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
> [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
Reply | Threaded
Open this post in threaded view
|

Re: 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 <
[hidden email]> 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 <[hidden email]>
> 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
>> [hidden email]
>> 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]]

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