Lettura shape file e visualizzazione layer e punti:

0. Librerie utilizzate maptools & rgdal
1. Con la funzione readOGR() leggo il file di formato .shp
2. Con le coordinate (che posso ricavare anche con la funzione coordinates()) posso aggiungere punti, aggiungere del testo oppure selezionare una parte del grafico
3. Per avere i dati relativa a questa parte occorre scaricare il file zip contenente i file al link e inserire al posto di PATH la posizione dove li salvate.

Esempio dati Austria

Dati point pattern - Esempio Lombardia

setwd("PATH")
require(maptools) 
require(rgdal)
# lettura shape file Lombardia in R (in coordinate GaussBoaga)
    lomb.poly<-readOGR("lombardia.shp",verbose=TRUE) 

OGR data source with driver: ESRI Shapefile 
Source: "PATH\lombardia.shp", layer: "lombardia"
with 1 features
It has 7 fields
Integer64 fields read as strings:  COM_ COM_ID 

    class(lomb.poly)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

    slotNames(lomb.poly)  # contenuto dell'oggetto

[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

# lettura dati 
    inc=incendi <- read.table("Incendi_2003.csv", header=T,sep=";");   
      str(incendi)        # descrizione sintetica

'data.frame':   380 obs. of  3 variables:
 $ Ettari: num  9 0.1 2 0.1 0.2 2 1.15 0.05 0.02 3 ...
 $ Nord  : int  5076900 5098880 5090650 5083300 5083300 5075000 5075200 5114790 5116640 4952750 ...
 $ EstN  : int  1604750 1518920 1508800 1491200 1491300 1563520 1638900 1601120 1603410 1524300 ...

# plot dati su mappa

    plot(lomb.poly)
    points(incendi$EstN,incendi$Nord)

Dati areali - Esempio Austria

setwd("PATH")
aus.poly<-readShapePoly("austrianuts3.shp", 
                        IDvar="ID",verbose=TRUE)            ### lettura shape file in R

Shapefile type: Polygon, (5), # of Shapes: 35

  slotNames(aus.poly)

[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

  slot(aus.poly,"data") # equivalente -> aus.poly@data  

         ID                     NAME
AT13   AT13                     Wien
AT111 AT111         Mittelburgenland
AT112 AT112           Nordburgenland
AT113 AT113            Südburgenland
AT121 AT121  Mostviertel-Eisenwurzen
AT122 AT122     Niederösterreich-Süd
AT123 AT123             Sankt Pölten
AT124 AT124              Waldviertel
AT125 AT125              Weinviertel
AT126 AT126   Wiener Umland/Nordteil
AT127 AT127    Wiener Umland/Südteil
AT211 AT211       Klagenfurt-Villach
AT212 AT212              Oberkärnten
AT213 AT213             Unterkärnten
AT221 AT221                     Graz
AT222 AT222                   Liezen
AT223 AT223  Östliche Obersteiermark
AT224 AT224            Oststeiermark
AT225 AT225  West- und Südsteiermark
AT226 AT226 Westliche Obersteiermark
AT311 AT311               Innviertel
AT312 AT312                Linz-Wels
AT313 AT313              Mühlviertel
AT314 AT314          Steyr-Kirchdorf
AT315 AT315             Traunviertel
AT321 AT321                   Lungau
AT322 AT322           Pinzgau-Pongau
AT323 AT323    Salzburg und Umgebung
AT331 AT331                Außerfern
AT332 AT332                Innsbruck
AT333 AT333                 Osttirol
AT334 AT334         Tiroler Oberland
AT335 AT335        Tiroler Unterland
AT341 AT341   Bludenz-Bregenzer Wald
AT342 AT342  Rheintal-Bodenseegebiet

  a=slot(aus.poly,"polygons"); 
  length(a); 

[1] 35

  str(a[[1]])

Formal class 'Polygons' [package "sp"] with 5 slots
  ..@ Polygons :List of 1
  .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. ..@ labpt  : num [1:2] 16.4 48.2
  .. .. .. ..@ area   : num 0.0502
  .. .. .. ..@ hole   : logi FALSE
  .. .. .. ..@ ringDir: int 1
  .. .. .. ..@ coords : num [1:335, 1:2] 16.2 16.2 16.2 16.2 16.2 ...
  ..@ plotOrder: int 1
  ..@ labpt    : num [1:2] 16.4 48.2
  ..@ ID       : chr "AT13"
  ..@ area     : num 0.0502

    plot(aus.poly,col = "green",border="red")   ### mappa regione 
    centr=coordinates(aus.poly)
    points(centr,col="red")
    text(centr[,1],centr[,2],aus.poly$ID,col="blue")
# estrazione sotto shape
    xx<-aus.poly[aus.poly$ID=="AT13",]          # poligono di Vienna
    plot(xx,add=T,col="blue")