עיבוד גיאוגרפי עם גאופאנדאס


גיאוגרפיה

הבסיס

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

ראשית, על מנת לפתוח קובץ נשתמש בפעולה,

import geopandas as gpd
police_stations_map = gpd.read_file('GeoLayers/PoliceStations/PoliceStations.shp')

הקובץ שנפתח יהיה מהסוג shape שכן לרוב שכבות שנוריד יגיעו בפורמט זה. אך, אדגיש כי הסיפרייה יודעת לעבוד גם עם סוגי קבצים אחרים כגון geojson,gpkg. הקובץ אותו פתחתי הוא קובץ של איזורי האחריות של כל תחנת משטרה, שנמצא באתר המאגרים הממשלתיים וזמין כאן, ובו אני אתעלל במדריך הזה.

ראשית, נראה שאם נבחן מה מכיל המשתנה אליו טענו את הקובץ,



נראה כי הוא בעקרון מתפקד כמו טבלת pandas רגילה (dataframe), עם שדה גיאומטרי. ושם המשתנה שלו הוא geodataframe. בעצם ניתן לבצע עליו את כל הפעולות שניתן על טבלה רגילה (query,apply,iloc) , בנוסף לפעולות הכוללות גיאומטריה, אותן נסקור כעת.

הצגת כמפה

הצגה פשוטה

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



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

  • figsize - גודל התמונה שתינתן, מתחיל ב(1,1). מצאתי כי הגודל (15,15) נותן תמונה ברזולוצייה טובה.
  • color - הצבע של האובייקטים במפה, מתקבל כסטרינג, לדוגמא 'red'.
  • edgecolors - צבע הגבול של פוליגון, מתקבל כסטרינג.
  • hatch - תבנית שתוכל על המפה, מבין התבניות הללו,

    התבנית מתקבלת כסטרינג, ואם נשים יותר מאותו הסימן, התבנית תהיה צפופה יותר. למשל '//' תיצור אלכסונים בצפיפות כפולה מהמוצג בתמונה.
  • linewidths - עובי גבול הפוליגון, מתקבל כמספר, מצאתי כי עובי סביר ניתן על ידי הספרה 1.
  • column - עמודה שאת המידה בה נרצה לייצג על המפה, לדוגמא - אם נרצה לייצג כמות אוכלוסייה בכל פוליגון, נבחר בעמודת כמות האוכלוסייה.
  • cmap - במידה ואופצייה זו בשימוש, הפוליגונים במפה ייצבעו לפי העמודה הנבחרת (column), לדוגמא, אם תבחר עמודת האוכלוסייה, הצבע הכי חזק יינתן לפוליגון בעל ערך האוכלוסייה הגבוהה ביותר. מיפוי הצבעים לפי חוזק ייבחר לפי הסטרינג שיוכנס המתאר מפת צבעים כלשהי מתוך המפות שנמצאות כאן, לדוגמא אני בחרתי במפה 'cool', כי ככה אנחנו, באתר הבקתה. אם נרצה לצבוע רנדומלית, כמו במפת העולם, ניתן לבחור את column כעמודת שם המקום, או משהו דומה, וכך הערך יחושב בצורה יחסית רנדומלית.
  • legend - האם נרצה להציג מקרא למפה, מקבל משתנה בוליאני.


וכעת התקבלה מפה שנראית הרבה יותר נורמלי.

הצגה כשכבות

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

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(15,15),facecolor='white')
police_stations_map.plot(ax = ax,cmap = 'gist_earth',zorder = 1,column = 'TahanaName', edgecolor = 'black', linewidths = 1)
schools_map.plot(ax = ax,zorder = 2,color = 'brown')

ראשית נייבא את הספרייה matplotlib ומתוכה את pyplot כplt. מהספרייה נשתמש בפעולה subplots שתחזיר לנו שני עצמים fig ו ax, הראשון מייצג את כלל התרשים, והשני מייצג את "מערכת הצירים" בתוך התרשים (שמכילה את במפה). ובעצם למשתנה השני, של מערכת הצירים, נוכל להוסיף מספר עצמים המכילים שכבות גיאוגרפיות ואז לראות אותן אחת על השנייה.
את ההוספה נבצע על ידי הגדרת המשתנה ax כשנבצע את פעולת הplot, המייצג את "מערכת הצירים" בה תשורטט המפה. אותו נגדיר כמערכת הצירים שיצרנו לפני, ששם המשתנה שלה הוא גם ax. ותתקבל המפה המאוחדת, לדוגמא כאן הוספתי שכבה של מוסדות חינוך על גבי המפה הקודמת.



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

מערכת ייחוס

למה?

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



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



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

מערכות ייחוס רלוונטיות בארץ

  • ITM - מערכת הייחוס הרשמית של ישראל היא ITM או בשמה המלא Israel Transverse Mercator המערכת נוצרה ב1993 והיא יחסית מודרנית, ובנויה על מנת לייצג את אזור הארץ בצורה מיטבית. הקאורדינטות בה נמדדות במטרים, כמרחק מנקודת האפס שנמצאת בערך באמצע הנילוס המצרי. כך שישראל בין 100,000 ל300,000 על ציר האורך, ובין 400,000 ל800,000 בציר הרוחב. לרוב המערכת נמצאת בשימוש במאגרים ממשלתיים. מזהה מערכת הייחוס שלה (EPSG) הוא 2039.
  • WSG84 - מערכת הייחוס הגלובאלית המקובלת ביותר, פותחה על ידי משרד הבטחון האמריקאי ב1984, ובנויה להציג את כלל כדור בארץ בצורה מיטבית (ככל שזה אפשרי). נמצאת בשימוש בגוגל מאפס, מערכת הלווין ורוב הגופים הבינלאומיים. הקאורדינטות בה מיוצגות על ידי קווי רוחב ואורך, כך שקווי הרוחב בין מינוס 90 בקוטב הדרומי ופלוס בצפוני, וקווי האורך בין מינוס 180 ממערב לארה"ב ופלוס ממזרח ליפן. ישראל נמצאת בה בין קו האורך 34 ל36, ועל קו הרוחב בין 29 ל33. מזהה מערכת הייחוס שלה (EPSG) הוא 4326.

נוכל לקבל מידע על מערכת ייחוס של קובץ שטענו בעזרת גישה לרכיב crs,



שייתן לנו מידע מפורט עליה, בנוסף למזהה שלה, הEPSG. בו נשתמש במידה ונרצה לבצע המרה למערכת ייחוס אחרת, למשל על מנת להמיר את שכבת אזורי האחריות של תחנות המשטרה לWGS84 מITM נבצע,

police_stations_map = police_stations_map.to_crs(4326)

טעינה מקובץ גולמי

במידה ונטען קובץ של אקסל או csv המערכת לא תדע אוטומטית איזו עמודה מתארת גיאומטריה, ובטח באיזו מערכת ייחוס הקובץ רשום, ונדרש לעשות כמה פעולות ידניות.
ראשית, נטען את הקובץ עם read_file, ומה שיתקבל יהיה geodataframe שעמודת הgeometry שלו ריקה. ואנו נדרש להשלים אותה בכל שורה לפי עמודות שמהן נגזור את המיקום. ישנן שתי אופציות מרכזיות לכך,

  • אם קיימות עמודות x,y בטבלה, מה שנעשה יהיה ראשית להבין באיזו מערכת ייחוס הן רשומות, אם המספרים באלפים ככה"נ זאת ITM הנמדדת במטרים, ואם קטנים ממאה זאת ככה"נ WSG84 הנמדדת לפי קווי אורך ורוחב. לאחר מכן נשתמש בפעולה הבאה,
    schools_map['geometry'] = gpd.points_from_xy(schools_map['X'], schools_map['Y'], crs="EPSG:2039")

    הפעולה תרשום בעמודה geometry את ערך הנקודה האחוד בהתבסס על שיעורים מהעמודות x,y, במערכת הייחוס "EPSG:2039" שהינה קוד ל ITM.
  • אם ישנה עמודת WKT (כלומר שמלאה בערכי Point או Polygon) נפתח את הקובץ כך,
    schools_map = gpd.read_file('schools_map.csv')
    schools_map['geometry'] = gpd.GeoSeries.from_wkt(schools_map['WKT'])
    schools_map = schools_map.set_crs(2039)

    בעצם נשתמש בפעולה gpd.GeoSeries.from_wkt() על מנת להמיר את ערכי הWKT לגיאומטריה, אותה נאחסן בעמודה geometry, המשומשת בפעולות הגיאומטריות על הטבלה. לאחר מכן, נזין את מערכת הייחוס המתאימה בעזרת set_crs(), והשכבה מוכנה לשימוש.

עיבוד גיאוגרפי

חישוב כיסוי

פעולת חישוב הכיסוי, overlay(), תחזיר לנו טבלה חדשה כתוצאה מחישוב כיסוי שנגדיר בין שתי שכבות. סוגי הכיסוי הינם,



על מנת להבין את אופי השימוש, בואו נקח דוגמא. אנחנו רוצים לשדוד בנק הנמצא במיקום לא ידוע בגוש 100270/1, ואנו רוצים להעריך מאיזו תחנת משטרה יגיעו השוטרים הבאים לנטרל אותנו. אך זה לא כזה פשוט שכן ארבעה אזורי אחריות שונים של תחנות משטרה נמצאים באותה הגוש. מה נעשה? נבדוק איזה אזור אחריות מכסה את מירב הגוש. נכתוב את הסקריפט הבא,

target_gush = gush_map.loc[(gush_map['GUSH_NUM'] == 100270) & (gush_map['GUSH_SUFFI'] == 1)]
police_in_gush = target_gush.overlay(police_stations_map,how='intersection')
police_in_gush['Area'] = police_in_gush.area
police_in_gush.plot(figsize=(15,15), edgecolor = 'black', linewidths = 1,column = 'Area',cmap = 'Reds')

ראשית, נאתר את פוליגון הגוש ונשמור אותו במשתנה target_gush לאחר מכן, ניקח את אותו הגוש ונבצע עליו פעולת overlay() עם מפת מחוזות המשטרה, במשתנה how נזין את סוג הכיסוי, ומכיוון שאנו רוצים את תחנות המשטרה שבתוך הגוש נבחר בintersection. המשתנה שיוחזר כעת יהיה פוליגונים של אותם אזורי האחריות, אך כל מה שמחוץ לגוש נחתך. לכן כעת נוכל לחשב את השטח של כל אזור אחריות, מה שנאחסן בעמודה Area. והתוצאה הינה,



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

חיבורים מרחביים

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



והאמת, שעל מנת לעשות זאת נשתמש רק בשורת קוד אחת,

schools_in_ashkelon = schools_in_ashkelon.sjoin(gushim_in_ashkelon, how="left")

נשתמש בפעולה sjoin (Spatial Join) על שכבת בתי הספר, ונזין אליה את שכבת הגושים, ועוד פרמטר שניגע בו תכף. מה שיתבצע כעת, זה שתיבדק כל נקודה בשכבת בתי הספר בנוגע לגוש בו היא נמצאת. והמידע המתאר את אותו הגוש משכבת הגושים יתווסף אליה. והטבלה (חלק ממנה) תראה כך,

NAME UNIQ_ID geometry GUSH_NUM GUSH_SUFFI SUB_GUSH_I
בית ספר - ממ"ד תורני ניצני קטיף 54449201 POINT (165309.320727 627014.844094) 2784.0 0.0 11401.0
ישיבה - בני יששכר 65109751 POINT (162614.418627 623124.147394) 3153.0 0.0 27446.0
בית ספר - מדעים 54446771 POINT (159910.121527 619896.065794) 2508.0 0.0 25646.0
גן ילדים - גן ממ"ד ניצן 54449973 POINT (158761.132827 618269.030594) 1965.0 0.0 11944.0
בית ספר - אור החיים 54451301 POINT (158256.037727 618153.504494) 1996.0 0.0 11949.0
בית ספר - תורני אסף מימון 54448980 POINT (161354.814027 619690.586594) 2383.0 0.0 19885.0
ישיבה - ע"ש הבן איש חי 54448028 POINT (158767.292427 618911.438894) 1962.0 0.0 11894.0
המרכז להכשרה טכנולוגית 54452647 POINT (161048.478227 618276.770994) 1207.0 0.0 11937.0
בית ספר - מוריה 54450273 POINT (158897.767327 618978.676294) 1963.0 0.0 11898.0
בית ספר - מ"מ אומנויות 54447000 POINT (160968.262227 619076.277994)
חווה לחינוך חקלאי אורות הסביבה ע"ש עלי בן צבי 54446961 POINT (158531.846727 619909.624294) 1942.0 0.0 11835.0

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

  • 'left' - הפעולה תקח את השורות מהטבלה עליה הופעלה הפעולה, ותמצא בתוך איזה פוליגון הן נמצאות מהטבלה השנייה. אם לא נמצאת התאמה, עדיין תשמר השורה.
  • 'right' - אותו הדבר כמו הפעולה הקודמת, פרט להיפוך בין הטבלאות. כלומר ישמרו כל השורות מהטבלה שהוזנה לפעולה כפרמטר.
  • 'inner' - הפעולה תשמור משתי הטבלאות רק שורות שיש ביניהן הצלבה.

חיבורים מרחביים ע"ב קירבה

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

gushim_in_ashkelon = gushim_in_ashkelon.sjoin_nearest(schools_in_ashkelon, how="left",distance_col='distance')
gushim_in_ashkelon.plot(cmap='Reds',column='distance',figsize=(15,15),legend=True)

הפעולה, sjoin_nearest(), מאוד דומה לחיבור הרגיל. הפרמטר הראשון יהיה השכבה שתוצלב, והפרמטר how עובד בצורה זהה. ולפרמטר distance_col נזין את שם העמודה בה ירשם המרחק הקצר ביותר לבית ספר. ניצור מהשכבה מפת חום והתוצאה תהיה,



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

best_sum = -1; best_schools_in_ashkelon = None
for index, row in gushim_in_ashkelon.iterrows():
____cur_schools_in_ashkelon = gpd.GeoDataFrame(pd.concat((schools_in_ashkelon, gpd.GeoDataFrame({'NAME':[row['GUSH_NUM']],'geometry':[row['geometry']]}, crs="EPSG:2039")), axis = 0))
____cur_gushim_in_ashkelon_join = gushim_in_ashkelon.sjoin_nearest(cur_schools_in_ashkelon, how="left",distance_col='distance')
____current_sum = cur_gushim_in_ashkelon_join.loc[cur_gushim_in_ashkelon_join['distance'] < 500].shape[0]
____if(current_sum > best_sum):
________best_sum = current_sum; best_schools_in_ashkelon = cur_schools_in_ashkelon

מה בעצם עשינו פה? אז ראשית הגדרנו שתי מתשתנים,

  • best_sum - יכיל את מספר הגושים הגדול ביותר שהגענו אליו, הנמצאים במרחק חצי קילומטר מבית ספר.
  • best_schools_in_ashkelon - יכיל את השכבה עם מיקום בית הספר האופטימלי.

וכעת נעבור בלולאה על כל אחד מהגושים באשקלון (gushim_in_ashkelon) , והלולאה תבצע את הפעולות הבאות (מסודר לפי שורה),

  1. ניצור משתנה זמני לאחסון שכבה שתכיל את בתי הספר, והגוש שהנבדק.
  2. נבצע חישוב מרחקים בין שכבת הגושים, והשכבה שיצרנו בסעיף הקודם.
  3. נרשום במשתנה את כמות הגושים שמרחקם קטן מחצי קילומטר עבור השכבה שיצרנו.
  4. במידה והכמות הכי גדולה עד עכשיו,
  5. הכמות, ושכבת בתי הספר עם המיקום החדש האידיאלי, יישמרו במשתנים המתאימים.

והתוצאה תהיה שבמשתנה best_schools_in_ashkelon יינתן מיקום בתי הספר בנוסף למיקום החדש האידיאלי, כך (המיקום החדש לא בצורת עיגול אלא כגוש),



עוד פרמטר שימושי שניתן להזין לפעולה sjoin_nearest() הוא max_distance שאליו נזין את המרחק המקסימלי אליו נסתכל בחיפוש אובייקט קרוב, כך שמעבר אליו החיפוש כבר אינו רלוונטי לשאילתה. יכלנו לדוגמא בעזרתו לחסוך כוח מחשוב באלגוריתם שרשמנו על ידי הגדרת המרחק כחצי קילומטר.

עוד מניפולציות

  • Dissolve - הפעולה משתמשת לאיחוד פוליגונים ע"ב פרמטר מסויים. לדוגמא, בשכבת אזורי אחריות תחנות המשטרה ישנה עמודה המשייכת את התחנה למחוז כלשהו, צפון, דרום וכו'. אם נרצה לגזור מהשכבה את השכבה של המחוזות נוכל להשתמש בפעולה,
    police_stations_map.dissolve(by="MahozName")

    אם לא נזין את הפרמטר by נקבל פשוט איחוד של כלל הקבוצה.
  • Centroid - הפעולה מקבלת גיאומטריות שונות ומחזירה את נקודת האמצע של כל אחת מהן. למשל אם נרצה שכבה של כל נקודות האמצע של אזורי האחריות,
    police_stations_map['geometry'] = police_stations_map['geometry'].centroid

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


אין תגובות:

הוסף רשומת תגובה

הירשם ל- רשומות (Atom)