כיצד לבצע ניתוח מרחבי ב- R עם sf

איפה אתה מצביע? מי אתם המחוקקים? מה המיקוד שלך? לשאלות אלו יש משהו במשותף מבחינה גיאו-מרחבית: התשובה כוללת קביעת איזה מצולע נקודה נופלת.

חישובים כאלה נעשים לעתים קרובות עם תוכנת GIS מיוחדת. אבל קל לעשות זאת גם ב- R. אתה צריך שלושה דברים:

  1. דרך לקידוד כתובות גיאוגרפיות למציאת קו רוחב ואורך; 
  2. קבצי צורה המתארים גבולות מצולע מיקוד; ו 
  3. חבילת ה- sf.

לצורך קידוד גיאוגרפי, אני בדרך כלל משתמש בממשק ה- API של geocod.io. זה בחינם ל -2,500 חיפושים ביום ויש לו חבילת R נחמדה, אבל אתה צריך מפתח API (בחינם) כדי להשתמש בו. כדי לעקוף את המורכבות הזו עבור מאמר זה, אשתמש בממשק ה- API של Nominatim עם קוד פתוח בחינם עם קוד פתוח. זה לא דורש מפתח. לחבילה tmaptools יש פונקציה,, geocode_OSM()להשתמש ב- API זה.

יבוא והכנה של נתונים גיאו-מרחביים

אני אשתמש בחבילות sf, tmaptools, tmap ו- dplyr. אם אתה רוצה לעקוב אחריו, טען כל אחד מהם עם pacman::p_load()או התקן אותו עדיין לא במערכת שלך install.packages(), ואז טען את כל אחד מהם library().

לדוגמא זו, אצור וקטור עם שתי כתובות, משרדנו בפרמינגהם, מסצ'וסטס, ומשרד RStudio בבוסטון.

כתובות <- c ("492 Path Connecticut Old, Framingham, MA",

"250 Northern Ave., בוסטון, MA")

קידוד גיאוגרפי הוא פשוט עם geocode_OSM. אתה יכול לראות את התוצאות על ידי הדפסת שלוש העמודות הראשונות כולל קו רוחב ואורך:

כתובות geocoded_ - geocode_OSM (כתובות)

הדפס (כתובות גיאוגרפיות [, 1: 3])

שאילתה lat lon

מס '1 492 נתיב קונטיקט הישן, פרמינגהם, MA 42.31348 -71.39105

# 2 250 Ave Ave, בוסטון, MA 42.34806 -71.03673

ישנן מספר דרכים להשיג קובצי צורת מיקוד. הקלה ביותר היא ככל הנראה אזורי טבלאות המיקוד של לשכת המפקד האמריקני, הדומים לגבולות שירות הדואר האמריקני אם לא ממש זהים.

ניתן להוריד קובץ ZCTA ישירות מלשכת האוכלוסין האמריקאית, אך זהו קובץ לכל הארץ. עשה זאת רק אם לא אכפת לך מקובץ נתונים גדול. 

מקום אחד להוריד קובץ ZCTA למדינה אחת הוא Census Reporter. חפש נתונים לפי מדינה, כגון אוכלוסייה, ואז הוסף מיקוד לגאוגרפיה ובחר הורדת נתונים כטופס צורה.

אני יכול לפתוח את הקובץ שהורדתי באופן ידני, אבל זה קל יותר ב- R. כאן אני משתמש unzip()בפונקציה של בסיס R בקובץ שהורד, ופותח אותו אל ספריית משנה לפרויקט בשם ma_zip_shapefile. junkpaths = TRUEהטיעון הזה אומר שאני לא רוצה לפתוח את רוכסן הוספת ספריית משנה נוספת המבוססת על שם קובץ ה- zip.

unzip ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", זבל מסלולים = TRUE,

להחליף = TRUE)

ייבוא ​​וניתוח גיאו-מרחבי עם sf

עכשיו סוף סוף עבודה גיאו-מרחבית. אני אייבא את הטופס לצורת R באמצעות st_read()הפונקציה של sf .

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # קריאת שכבה 'acs2017_5yr_B01003_86000US02648' ממקור הנתונים `/ משתמשים /smachlis/Documents/MoreWithR/ma_zip_shp_shp תכונות וארבעה שדות # סוג גיאומטריה: MULTIPOLYGON # ממד: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs

צירפתי את תגובת המסוף בזמן ההפעלה st_read()מכיוון שמוצג שם מידע כלשהו: epsg. זה אומר באיזו מערכת הפניה לקואורדינטות שימש ליצירת הקובץ . כאן זה היה 4326. מבלי להיכנס עמוק מדי לעשבים שוטפים, epsg בעצם מציין  איזו מערכת שימשה לתרגום אזורים על גבי כדור תלת מימד - כדור הארץ - לקואורדינטות דו ממדיות (קו רוחב ואורך) . זה חשוב מכיוון שיש הרבה מערכות התייחסות שונות לקואורדינטות. אני רוצה שפוליגוני המיקוד ונקודות הכתובת שלי ישתמשו באותו, כך שהם יתייצבו כראוי.

הערה: קובץ זה כולל במקרה מצולע לכל מדינת מסצ'וסטס, שאיני זקוק לו. אז אני אסנן את השורה הזו במסצ'וסטס עם

zipcode_geo <- dplyr :: filter (מיקוד_geo,

שם! = "מסצ'וסטס")

מיפוי טופס הצורה בעזרת ממפה

מיפוי נתוני המצולע אינו הכרחי, אך זה בדיקה יפה של קובץ הצורה שלי כדי לראות אם הגיאומטריה היא מה שאני מצפה. אתה יכול לבצע עלילה מהירה של אובייקט sf עם הפונקציה tmap qtm()(קיצור של מפת נושאים מהירה).

qtm (zipcode_geo) +

tm_legend (show = FALSE)

מסכים שצילם שרון מכליס,

וזה נראה כאילו אכן יש לי גיאומטריה של מסצ'וסטס עם מצולעים שיכולים להיות מיקודים.

בהמשך אני רוצה להשתמש בנתוני הכתובת הגיאוגרפית. זוהי כרגע מסגרת נתונים רגילה, אך יש להמיר אותה לאובייקט geospatial sf עם מערכת הקואורדינטות הנכונה.

אנו יכולים לעשות זאת בעזרת st_as_sf()הפונקציה של sf . (הערה: פונקציות חבילה sf הפועלות על נתונים מרחביים מתחילות עם st_, מה שמייצג "מרחבי" ו"זמני ".)

st_as_sf()לוקח כמה ויכוחים. בקוד שלמטה, הטיעון הראשון הוא האובייקט להמרה - הכתובות המקודדות שלי. וקטור הארגומנט השני אומר לפונקציה באילו עמודות יש את הערכים x (קו אורך) ו- y (קו רוחב). השלישי מגדיר את מערכת ההפניה לקואורדינטות ל- 4326, כך שהיא זהה לפוליגוני המיקוד שלי.

point_geo <- st_as_sf (כתובות GeoCoded_,

קואורדיות = c (x = "lon", y = "lat"),

crs = 4326)

Geospatial מצטרף עם sf

עכשיו, כשהגדרתי את שתי ערכות הנתונים שלי, חישוב המיקוד עבור כל כתובת קל בעזרת st_join()הפונקציה של sf . התחביר:

st_join (point_sf_object, polygon_sf_object, join = join_type)

בדוגמה זו, אני רוצה לרוץ st_join()על הנקודות המקודדות גיאוגרפית ראשון ועל מצולעי המיקוד השני. זה מה שנקרא פורמט הצטרפות שמאלי: כל הנקודות בנתונים הראשונים (כתובות מקודדות גיאוגרפיות) כלולות, אך רק נקודות בנתונים השנייה (מיקוד) שתואמות. לבסוף, סוג ההצטרפות שלי הוא st_withinמכיוון שאני רוצה שההתאמה תהיה נקודות בתוכה. 

התוצאות שלי <- st_join (point_geo, מיקוד_geo,

הצטרף = st_within)

זהו זה! עכשיו אם אני מסתכל על התוצאות שלי על ידי הדפסת כמה מהעמודות החשובות ביותר, תראה שלכל כתובת יש מיקוד (בעמודה "שם"). 

הדפס (my_results [, c ("query", "name", "geometry")])

# אוסף תכונות פשוט עם 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 Path Old Connecticut, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2 250 Northern Ave., בוסטון, MA 02210 POINT (-71.03673 42.34806)

מיפוי נקודות ומצולעים עם ממפה

אם תרצה למפות את הנקודות והמצולעים, הנה דרך אחת לעשות זאת באמצעות tmap:

tm_shape (מיקוד_geo) +

tm_fill () +

tm_shape (התוצאות שלי) +

tm_bubbles (col = "red", size = 0.25)

צילום מסך של שרון מכליס,

רוצה עוד טיפים ל- R? עבור לדף "עשה יותר עם R"!