Πώς να κάνετε χωρική ανάλυση στο R με sf

Πού ψηφίζετε; Ποιοι είστε νομοθέτες; Ποιος είναι ο ταχυδρομικός κώδικας σας; Αυτές οι ερωτήσεις έχουν κάτι κοινό από γεωγραφική άποψη: Η απάντηση περιλαμβάνει τον προσδιορισμό του πολυγώνου σε ένα σημείο.

Τέτοιοι υπολογισμοί γίνονται συχνά με εξειδικευμένο λογισμικό GIS. Αλλά είναι επίσης εύκολο να το κάνετε στο R. Χρειάζεστε τρία πράγματα:

1. Ένας τρόπος για να διευθύνσεις γεωκωδικού για εύρεση γεωγραφικού πλάτους και μήκους.
2. Shapefiles που περιγράφουν τα όρια πολυγώνων ταχυδρομικού κώδικα. και
3. Το πακέτο sf.

Για τη γεωκωδικοποίηση, συνήθως χρησιμοποιώ το API geocod.io. Είναι δωρεάν για 2.500 αναζητήσεις την ημέρα και έχει ένα ωραίο πακέτο R, αλλά χρειάζεστε ένα (δωρεάν) κλειδί API για να το χρησιμοποιήσετε. Για να ξεπεράσω αυτό το κομμάτι της πολυπλοκότητας αυτού του άρθρου, θα χρησιμοποιήσω το δωρεάν, ανοιχτού κώδικα Open Street Map Nominatim API. Δεν απαιτεί κλειδί. Το πακέτο tmaptools έχει μια λειτουργία, geocode_OSM (), για να χρησιμοποιήσετε αυτό το API.

Εισαγωγή και προετοιμασία γεωχωρικών δεδομένων

Θα χρησιμοποιώ τα πακέτα sf, tmaptools, tmap και dplyr. Εάν θέλετε να ακολουθήσετε, φορτώστε το κάθε ένα με pacman :: p_load () ή εγκαταστήστε κανένα ακόμη στο σύστημά σας με install.packages ()και στη συνέχεια φορτώστε το κάθε ένα με βιβλιοθήκη().

Για αυτό το παράδειγμα, θα δημιουργήσω έναν φορέα με δύο διευθύνσεις, το γραφείο μας στο Framingham της Μασαχουσέτης και το γραφείο RStudio στη Βοστώνη.

διευθύνσεις <- c ("492 Old Connecticut Path, Framingham, MA",
"250 Northern Ave., Boston, MA")

Η γεωκωδικοποίηση είναι απλή με το geocode_OSM. Μπορείτε να δείτε τα αποτελέσματα εκτυπώνοντας τις τρεις πρώτες στήλες συμπεριλαμβανομένου του γεωγραφικού πλάτους και μήκους:

geocoded_addresses <- geocode_OSM (διευθύνσεις)
εκτύπωση (geocoded_addresses [, 1: 3])
ερώτημα lat lon
# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105
# 2 250 Northern Ave., Βοστόνη, MA 42.34806 -71.03673

Υπάρχουν διάφοροι τρόποι για να λάβετε shapefiles ταχυδρομικού κώδικα. Το πιο εύκολο είναι πιθανότατα οι Περιοχές Πίνακας Ταχυδρομικού Κώδικα του Γραφείου Απογραφής των ΗΠΑ, οι οποίες μοιάζουν με, αν όχι ακριβώς οι ίδιες με τα όρια των ταχυδρομικών υπηρεσιών των ΗΠΑ.

Μπορείτε να κατεβάσετε ένα αρχείο ZCTA απευθείας από το Γραφείο Απογραφής των ΗΠΑ, αλλά είναι ένα αρχείο για ολόκληρη τη χώρα. Κάνετε αυτό μόνο εάν δεν σας πειράζει ένα μεγάλο αρχείο δεδομένων.

Ένα μέρος για τη λήψη ενός αρχείου ZCTA για μία κατάσταση είναι το Census Reporter. Αναζητήστε οποιαδήποτε δεδομένα κατά πολιτεία, όπως πληθυσμό και, στη συνέχεια, προσθέστε ταχυδρομικό κώδικα στη γεωγραφία και επιλέξτε λήψη δεδομένων ως shapefile.

Θα μπορούσα να αποσυμπιέσω το ληφθέν αρχείο με μη αυτόματο τρόπο, αλλά είναι πιο εύκολο στο R. Εδώ χρησιμοποιώ το βασικό R ανοίγω φερμουάρ() λειτουργία σε ένα ληφθέν αρχείο και αποσυμπιέστε το σε έναν υποκατάλογο έργου που ονομάζεται ma_zip_shapefile. Οτι junkpaths = ΑΛΗΘΟΣ Το επιχείρημα λέει ότι δεν θέλω να αποσυμπιέσω προσθέτοντας έναν άλλο υποκατάλογο με βάση το όνομα του αρχείου zip.

αποσυμπιέστε ("data / acs2017_5yr_B01003_86000US02648.zip",
exdir = "ma_zip_shapefile", junkpaths = TRUE,
αντικατάσταση = ΑΛΗΘΟΣ)

Γεωχωρική εισαγωγή και ανάλυση με sf

Τώρα επιτέλους κάποια γεωχωρική εργασία. Θα εισαγάγω το shapefile στο R χρησιμοποιώντας sf's st_read () λειτουργία.

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Επίπεδο ανάγνωσης `acs2017_5yr_B01003_86000US02648 'από την πηγή δεδομένων` χαρακτηριστικά και 4 πεδία # τύπος γεωμετρίας: MULTIPOLYGON # διάσταση: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_d

Έχω συμπεριλάβει την απόκριση της κονσόλας κατά την εκτέλεση st_read () επειδή υπάρχουν κάποιες πληροφορίες που εμφανίζονται εκεί: το epsg. Που λέει ποιο σύστημα αναφοράς συντεταγμένων χρησιμοποιήθηκε για τη δημιουργία του αρχείου. Εδώ ήταν 4326. Χωρίς να μπαίνουμε πολύ βαθιά στα ζιζάνια, ένα βασικό μήνυμα δείχνειποιο σύστημα χρησιμοποιήθηκε για τη μετάφραση περιοχών σε έναν τρισδιάστατο κόσμο - τη Γη - σε δισδιάστατες συντεταγμένες (γεωγραφικό πλάτος και μήκος). Αυτό είναι σημαντικό επειδή υπάρχουν παρτίδα διαφορετικών συστημάτων αναφοράς συντεταγμένων. Θέλω τα πολύγωνα και τα σημεία διεύθυνσής μου να χρησιμοποιούν το ίδιο, ώστε να ευθυγραμμίζονται σωστά.

Σημείωση: Αυτό το αρχείο περιλαμβάνει ένα πολύγωνο για ολόκληρη την πολιτεία της Μασαχουσέτης, το οποίο δεν χρειάζομαι. Θα φιλτράρω λοιπόν αυτήν τη σειρά της Μασαχουσέτης

zipcode_geo <- dplyr :: φίλτρο (zipcode_geo,
όνομα! = "Μασαχουσέτη")

Χαρτογράφηση του shapefile με tmap

Η χαρτογράφηση των δεδομένων πολυγώνου δεν είναι απαραίτητη, αλλά είναι ένας καλός έλεγχος του shapefile μου για να δω αν η γεωμετρία είναι αυτό που περιμένω. Μπορείτε να κάνετε μια γρήγορη πλοκή ενός αντικειμένου sf με tmap's qtm () (συντόμευση για γρήγορο χάρτη θέματος) λειτουργία.

qtm (zipcode_geo) +
tm_legend (εμφάνιση = FALSE)

Οθόνες που πυροβολήθηκαν από τον Sharon Machlis,

Και φαίνεται ότι έχω πράγματι γεωμετρία της Μασαχουσέτης με πολύγωνα που θα μπορούσαν να είναι ταχυδρομικοί κώδικες.

Στη συνέχεια θέλω να χρησιμοποιήσω τα δεδομένα διεύθυνσης με γεωκωδικοποίηση. Αυτή τη στιγμή είναι ένα απλό πλαίσιο δεδομένων, αλλά πρέπει να μετατραπεί σε ένα γεωγραφικό αντικείμενο sf με το σωστό σύστημα συντεταγμένων.

Μπορούμε να το κάνουμε με sf's st_as_sf () λειτουργία. (Σημείωση: ξεκινούν οι λειτουργίες πακέτου sf που λειτουργούν σε χωρικά δεδομένα st_, που σημαίνει «χωρικά» και «χρονικά».)

st_as_sf () παίρνει πολλά επιχειρήματα. Στον παρακάτω κώδικα, το πρώτο όρισμα είναι το αντικείμενο που πρέπει να μετασχηματιστεί - οι γεωκωδικοποιημένες διευθύνσεις μου. Το δεύτερο διάνυσμα ορίσματος λέει τη συνάρτηση ποιες στήλες έχουν τις τιμές x (γεωγραφικό μήκος) και y (γεωγραφικό πλάτος). Το τρίτο θέτει το σύστημα αναφοράς συντεταγμένων σε 4326, έτσι είναι το ίδιο με τα πολύγωνα ταχυδρομικού κώδικα.

point_geo <- st_as_sf (geocoded_addeses,
coords = c (x = "lon", y = "lat"),
crs = 4326)

Το Geospatial ενώνεται με το sf

Τώρα που έχω δημιουργήσει τα δύο σύνολα δεδομένων μου, ο υπολογισμός του ταχυδρομικού κώδικα για κάθε διεύθυνση είναι εύκολος με το sf st_join () λειτουργία. Η σύνταξη:

st_join (point_sf_object, polygon_sf_object, join = join_type)

Σε αυτό το παράδειγμα, θέλω να τρέξω st_join () πρώτα στα γεωκωδικοποιημένα σημεία και το δεύτερο ταχυδρομικό κώδικα πολυγώνων. Είναι μια λεγόμενη αριστερή ένωση: Ολα Περιλαμβάνονται σημεία στα πρώτα δεδομένα (διευθύνσεις με γεωκωδικοποίηση), αλλά μόνο σημεία στο δεύτερο (ταχυδρομικός κώδικας) δεδομένα που ταιριάζουν Τέλος, ο τύπος συμμετοχής μου είναι st_within, αφού θέλω ο αγώνας να είναι βαθμοί εντός.

my_results <- st_join (point_geo, zipcode_geo,
εγγραφή = st_within)

Αυτό είναι! Τώρα αν κοιτάξω τα αποτελέσματά μου εκτυπώνοντας αρκετές από τις πιο σημαντικές στήλες, θα δείτε ότι κάθε διεύθυνση έχει ταχυδρομικό κώδικα (στη στήλη "όνομα").

εκτύπωση (my_results [, c ("query", "name", "geometry")])
# Απλή συλλογή χαρακτηριστικών με 2 δυνατότητες και 2 πεδία # τύπος γεωμετρίας: ΣΗΜΕΙΟ # διάσταση: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # γεωμετρία ονόματος ερωτήματος # 1 492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2 250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Σημεία χαρτογράφησης και πολύγωνα με tmap

Αν θέλετε να χαρτογραφήσετε τα σημεία και τα πολύγωνα, ακολουθήστε έναν τρόπο για να το κάνετε με το tmap:

tm_shape (zipcode_geo) +
tm_fill () +
tm_shape (my_results) +
tm_bubbles (col = "κόκκινο", μέγεθος = 0,25)

Στιγμιότυπο από τον Sharon Machlis,

Θέλετε περισσότερες συμβουλές R; Μεταβείτε στη σελίδα "Κάντε περισσότερα με R"!