Comment faire une analyse spatiale en R avec sf

Où votez-vous? Qui êtes-vous les législateurs? Quel est votre code postal? Ces questions ont quelque chose de géospatial en commun: la réponse consiste à déterminer dans quel polygone se trouve un point.

Ces calculs sont souvent effectués avec un logiciel SIG spécialisé. Mais ils sont également faciles à faire dans R. Vous avez besoin de trois choses:

  1. Un moyen de géocoder les adresses pour trouver la latitude et la longitude; 
  2. Fichiers de formes qui délimitent les limites des polygones de code postal et 
  3. Le paquet sf.

Pour le géocodage, j'utilise généralement l'API geocod.io. Il est gratuit pour 2500 recherches par jour et dispose d'un joli package R, mais vous avez besoin d'une clé API (gratuite) pour l'utiliser. Pour contourner ce peu de complexité de cet article, j'utiliserai l'API Open Street Map Nominatim gratuite et open source. Il ne nécessite pas de clé. Le package tmaptools a une fonction,, geocode_OSM()pour utiliser cette API.

Importation et préparation de données géospatiales

J'utiliserai les packages sf, tmaptools, tmap et dplyr. Si vous voulez suivre, chargez chacun avec pacman::p_load()ou installez-en tout pas encore sur votre système avec install.packages(), puis chargez chacun avec library().

Pour cet exemple, je vais créer un vecteur avec deux adresses, notre bureau à Framingham, Massachusetts et le bureau RStudio à Boston.

adresses <- c ("492 Old Connecticut Path, Framingham, MA",

«250 Northern Ave., Boston, MA»)

Le géocodage est simple avec geocode_OSM. Vous pouvez voir les résultats en imprimant les trois premières colonnes, y compris la latitude et la longitude:

geocoded_addresses <- geocode_OSM (adresses)

print (adresses_ géocodées [, 1: 3])

requête lat lon

# 1492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105

# 2250 Northern Ave., Boston, MA 42.34806 -71.03673

Il existe plusieurs façons d'obtenir des fichiers de formes de code postal. Le plus simple est probablement les zones de tabulation des codes postaux du US Census Bureau, qui sont similaires sinon exactement les mêmes que les limites du US Postal Service.

Vous pouvez télécharger un fichier ZCTA directement à partir du US Census Bureau, mais c'est un fichier pour tout le pays. Ne faites cela que si cela ne vous dérange pas un gros fichier de données. 

Un endroit pour télécharger un fichier ZCTA pour un seul État est Census Reporter. Recherchez les données par état, telles que la population, puis ajoutez le code postal à la géographie et choisissez télécharger les données sous forme de fichier de formes.

Je pourrais décompresser mon fichier téléchargé manuellement, mais c'est plus facile dans R. Ici, j'utilise la fonction de base R unzip()sur un fichier téléchargé, et je le décompresse dans un sous-répertoire de projet nommé ma_zip_shapefile. Cet junkpaths = TRUEargument dit que je ne veux pas décompresser en ajoutant un autre sous-répertoire basé sur le nom du fichier zip.

unzip ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpaths = TRUE,

écraser = VRAI)

Importation et analyse géospatiales avec SF

Maintenant enfin quelques travaux géospatiaux. Je vais importer le fichier de formes dans R en utilisant la st_read()fonction de sf .

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Lecture de la couche `acs2017_5yr_B01003_86000US02648 'à partir de la source de données` /Users/smachlis/Documents_B01003_86000US02648. caractéristiques et 4 champs # type de géométrie: MULTIPOLYGON # dimension: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs

J'ai inclus la réponse de la console lors de l'exécution st_read()car certaines informations y sont affichées: le fichier epsg. Cela indique quel système de référence de coordonnées a été utilisé pour créer le fichier . Ici, c'était 4326. Sans aller trop loin dans les mauvaises herbes, un epsg indique essentiellement  quel système a été utilisé pour traduire des zones d'un globe tridimensionnel - la Terre - en coordonnées bidimensionnelles (latitude et longitude) . Ceci est important car il existe de nombreux systèmes de référence de coordonnées différents. Je veux que mes polygones de code postal et mes points d'adresse utilisent le même, afin qu'ils s'alignent correctement.

Remarque: ce fichier inclut un polygone pour tout l'état du Massachusetts, dont je n'ai pas besoin. Alors je vais filtrer cette ligne du Massachusetts avec

zipcode_geo <- dplyr :: filter (zipcode_geo,

name! = "Massachusetts")

Mappage du fichier de formes avec tmap

La cartographie des données de polygones n'est pas nécessaire, mais c'est une belle vérification de mon fichier de formes pour voir si la géométrie est ce que j'attends. Vous pouvez faire un tracé rapide d'un objet sf avec la fonction tmap qtm()(abréviation de carte de thème rapide).

qtm (zipcode_geo) +

tm_legend (afficher = FALSE)

Écrans tournés par Sharon Machlis,

Et il semble que j'ai effectivement une géométrie du Massachusetts avec des polygones qui pourraient être des codes postaux.

Ensuite, je veux utiliser les données d'adresse géocodées. Il s'agit actuellement d'un bloc de données simple, mais il doit être converti en un objet géospatial sf avec le bon système de coordonnées.

Nous pouvons le faire avec la st_as_sf()fonction de sf . (Remarque: les fonctions du package sf qui fonctionnent sur des données spatiales commencent par st_, qui signifie «spatial» et «temporel».)

st_as_sf()prend plusieurs arguments. Dans le code ci-dessous, le premier argument est l'objet à transformer - mes adresses géocodées. Le deuxième vecteur d'argument indique à la fonction quelles colonnes ont les valeurs x (longitude) et y (latitude). Le troisième définit le système de référence de coordonnées sur 4326, donc c'est le même que mes polygones de code postal.

point_geo <- st_as_sf (adresses_ géocodées,

coords = c (x = "lon", y = "lat"),

crs = 4326)

Geospatial s'associe à SF

Maintenant que j'ai configuré mes deux ensembles de données, le calcul du code postal pour chaque adresse est facile avec la st_join()fonction sf . La syntaxe:

st_join (point_sf_object, polygon_sf_object, join = join_type)

In this example, I want to run st_join() on the geocoded points first and the ZIP code polygons second. It’s a so-called left join format: All points in the first data (geocoded addresses) are included, but only points in the second (ZIP code) data that match. Finally, my join type is st_within, since I want the match to be points within. 

my_results <- st_join(point_geo, zipcode_geo,

join = st_within)

That’s it! Now if I look at my results by printing out several of the most important columns, you”ll see each address has a ZIP code (in the “name” column). 

print(my_results[,c("query", "name", "geometry")])

# Collection d'entités simples avec 2 caractéristiques et 2 champs # type de géométrie: POINT # dimension: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # nom de la requête géométrie # 1492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Cartographie des points et des polygones avec tmap

Si vous souhaitez mapper les points et les polygones, voici une façon de le faire avec tmap:

tm_shape (zipcode_geo) +

tm_fill () +

tm_shape (mes_résultats) +

tm_bubbles (col = "rouge", taille = 0,25)

Capture d'écran de Sharon Machlis,

Vous voulez plus de conseils R? Rendez-vous sur la page «Faites plus avec R»!