Finding Centroid point from shape file and labelling map



Finding Centroid point from shape file and labelling map


Finding Centroid point from shape file and map data labelling

In this test I try to find centroid point in shp file. In additionally I will show how to use this centroid point information when to labelling area name into map.
R Library PBSmapping offers wide range of functions to do that job.

Needed Library

require(rgdal)
## Loading required package: rgdal
## Loading required package: sp
## Warning: package 'sp' was built under R version 2.15.1
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.0, released 2011/12/29 Path to GDAL shared
## files: C:/Program Files/R/R-2.14.2/library/rgdal/gdal Loaded PROJ.4
## runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470] Path to PROJ.4
## shared files: C:/Program Files/R/R-2.14.2/library/rgdal/proj
require(PBSmapping)
## Loading required package: PBSmapping
## ----------------------------------------------------------- PBS Mapping
## 2.62.34 -- Copyright (C) 2003-2012 Fisheries and Oceans Canada
## 
## PBS Mapping comes with ABSOLUTELY NO WARRANTY; for details see the file
## COPYING. This is free software, and you are welcome to redistribute it
## under certain conditions, as outlined in the above file.
## 
## A complete user guide 'PBSmapping-UG.pdf' is located at C:/Program
## Files/R/R-2.14.2/library/PBSmapping/doc/PBSmapping-UG.pdf
## 
## To see demos, type '.PBSfigs()'.
## 
## Packaged on 2012-03-01 Pacific Biological Station, Nanaimo
## -----------------------------------------------------------
require(maptools)
## Loading required package: maptools
## Loading required package: foreign
## Loading required package: lattice
## Checking rgeos availability: TRUE

Importing Shapefile

FIN_adm2.pb <- importShapefile("FIN_adm2")

ImportShapefile reads the .prj file if it exists, but it
does not adopt the proj4 format

proj.abbr <- attr(FIN_adm2.pb, "projection")
# abbreviated projection info

print(proj.abbr)

Generate map using PBSmapping plotting functions

plotPolys(FIN_adm2.pb, projection = proj.abbr, 
    border = "gray", xlab = "Longitude", ylab = "Latitude")

plot of chunk unnamed-chunk-4

Calculate the centroids of adm2 area

FIN.centroids2 <- calcCentroid(FIN_adm2.pb, rollup = 1)

Check results

str(FIN.centroids2)
## Classes 'PolyData' and 'data.frame': 21 obs. of  3 variables:
##  $ PID: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ X  : num  29.9 27.5 26 27.9 26.3 ...
##  $ Y  : num  62.9 63.2 61.4 61.9 67.7 ...
##  - attr(*, "projection")= chr "LL"

Plot the centroids against the source polygons

In this purpose we use the PBSMapping package graphics routines. If you want you can plot these centroid point in to the map

plotPolys(FIN_adm2.pb, projection = proj.abbr, 
    border = "gray", xlab = "Longitude", ylab = "Latitude")
addPoints(FIN.centroids2, col = "blue", pch = 19)

plot of chunk unnamed-chunk-7

check the variables from imported shape file

str(FIN_adm2.pb)
## Classes 'PolySet' and 'data.frame':  136271 obs. of  5 variables:
##  $ PID: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ SID: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ POS: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ X  : num  28.5 28.5 28.5 28.5 28.6 ...
##  $ Y  : num  63.9 63.9 63.9 63.9 63.9 ...
##  - attr(*, "PolyData")=Classes 'PolyData' and 'data.frame':  21 obs. of  12 variables:
##   ..$ PID      : int  1 2 3 4 5 6 7 8 9 10 ...
##   ..$ ID_0     : int  78 78 78 78 78 78 78 78 78 78 ...
##   ..$ ISO      : Factor w/ 1 level "FIN": 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ NAME_0   : Factor w/ 1 level "Finland": 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ ID_1     : int  1 1 1 1 2 3 3 4 4 4 ...
##   ..$ NAME_1   : Factor w/ 5 levels "Eastern Finland",..: 1 1 1 1 2 3 3 4 4 4 ...
##   ..$ ID_2     : int  1 2 3 4 5 6 7 8 9 10 ...
##   ..$ NAME_2   : Factor w/ 19 levels "Central Finland",..: 8 9 13 17 7 5 10 3 6 13 ...
##   ..$ NL_NAME_2: Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
##   ..$ VARNAME_2: Factor w/ 19 levels "Etelä-Karjala|Södra Karelen",..: 12 14 15 3 10 5 13 4 9 15 ...
##   ..$ TYPE_2   : Factor w/ 1 level "Maakunta|landskap": 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ ENGTYPE_2: Factor w/ 1 level "Region": 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "parent.child")= num  1 1 1 1 0 1 0 0 0 0 ...
##  - attr(*, "shpType")= int 5
##  - attr(*, "prj")= __truncated__
##  - attr(*, "projection")= chr "LL"

To read dbf file we need foreign library

library(foreign)

reading shape file data from dbf file (ESRI files include normally: shp, dbf, prj, csv and shx data files)

labelData <- read.dbf("FIN_adm2.dbf", as.is = TRUE)

Data checking

head(labelData)
##   ID_0 ISO  NAME_0 ID_1          NAME_1 ID_2            NAME_2 NL_NAME_2
## 1   78 FIN Finland    1 Eastern Finland    1     North Karelia      <NA>
## 2   78 FIN Finland    1 Eastern Finland    2     North Savonia      <NA>
## 3   78 FIN Finland    1 Eastern Finland    3 Päijänne Tavastia      <NA>
## 4   78 FIN Finland    1 Eastern Finland    4  Southern Savonia      <NA>
## 5   78 FIN Finland    2         Lapland    5           Lapland      <NA>
## 6   78 FIN Finland    3            Oulu    6            Kainuu      <NA>
##                         VARNAME_2            TYPE_2 ENGTYPE_2
## 1   Pohjois-Karjala|Norra Karelen Maakunta|landskap    Region
## 2      Pohjois-Savo|Norra Savolax Maakunta|landskap    Region
## 3 Päijät-Häme|Päijänne Tavastland Maakunta|landskap    Region
## 4        Etelä-Savo|Södra Savolax Maakunta|landskap    Region
## 5                  Lappi|Lappland Maakunta|landskap    Region
## 6               Kainuu|Kajanaland Maakunta|landskap    Region

merging shape files area name and centroid points (x and y coordinates) labelData as data frame

labelData <- as.data.frame(labelData)

merging labelData and centroid data table fields

label <- data.frame(labelData$ID_2, FIN.centroids2$X, 
    FIN.centroids2$Y, labelData$NAME_2)

Event data function needs right column names

names(label) <- c("EID", "X", "Y", "label")

Note: Above routine used in eventdata function

#######names(label)<- c(“PID”, “X”, “Y”, “label”) –> USED IN POLYDATA FUNCTION

do the right data format into labelling

data <- as.EventData(label, projection = NULL, 
    zone = NULL)

Labelling into map

#col = 1 –> font color 1=black, 2=….
#font = 3 –> font size
#cex = 0.5 –>expansion factor

plotPolys(FIN_adm2.pb, projection = proj.abbr, 
    border = "gray", xlab = "Longitude", ylab = "Latitude")

addLabels(data, placement = "DATA", polys = FIN_adm2.pb, 
    col = 1, font = 3, cex = 0.5)
## Warning: The data's 'projection' attribute (NULL) differs from the
## projection of the plot region (LL).

plot of chunk unnamed-chunk-16

Have fun with mapping and RStudio

Marko

ps. This post is produced by using RStudio and knitr. Thanks to Yihui Xie for amazing work with this issue.

Copyright MySci 2019
Tech Nerd theme designed by Siteturner