יותר

המרת קואורדינטות למסדר מפקד עבור מערך נתונים גדול

המרת קואורדינטות למסדר מפקד עבור מערך נתונים גדול


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

הסתכלתי על השרשור הזה: ממשק API חינם להחזרת קו רוחב גיאוגרפי, קו אורך למערכת המפקד? עם זאת, אני חושב שכל פתרון http / API יהיה איטי מדי. כמו כן, אני לא חושב ש- QGIS יכול להתמודד עם כל כך הרבה נקודות.

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


אם אתה רוצה שזה יהיה מהיר ויש לך מספר רב של lat / longs לבדוק אני ממליץ בחום להשתמש ב- postgresql / PostGIS. עם זאת, ניתן לעשות זאת בפייתון, בתנאי שהתקינו חבילות OGR / GDAL כהלכה.

להלן מתאר כיצד ניתן לעשות זאת אך ורק בפייתון, בהנחה שהתקנת כראוי חבילות OGR / GDAL (התקנות אלה יכולות להיות לא טריוויאליות, אז היזהר)

מ osgeo יבוא ogr # מניח שהנקודות והצורה שלך כבר באותו נתון / השלכה shapefile_name = "census_tracts.shp" # גרסה זו לוקחת רשימה ארוכה_למטה של ​​הטופס שלמטה ושם צורה קובץ def getCensusTracts (long_lat_list, shapefile_name): driver = ogr .GetDriverByName ("ESRI Shapefile") dataSource = driver.Open (shapefile_name, 0) layer = dataSource.GetLayer () results_dict = {} i = 0 for feature in layer: geom = feature.GetGeometryRef () i + = 1 for pt ב- long_lat_list: gid = pt [0] lon = pt [1] lat = pt [2] point = ogr. Geometry (ogr.wkbPoint) point.AddPoint (lon, lat) if point. בתוך (geom) == נכון: feat_id = feature.GetField ("fips") אם gid בתוצאות_דיקה ו- feat_id לא בתוצאות_דיקט [gid]: results_dict [gid] .append (feat_id) אחר: results_dict [gid] = [feat_id] עבור pt ב- long_lat_list: gid = pt [0] lon = pt [1] lat = pt [2] if gid not in results_dict: results_dict [gid] = ['NA'] return results_dict # היכן האלמנטים נמצאים [id, long, lat] long_lat_list = [[1, -87.5, 35.5], [2, -78. 446, 41.353]] results_dict = getCensusTracts (long_lat_list, shapefile_name) #results_dict מחזיר מילון שבו {'id: list_of_fips_codes} מדפיס תוצאות_דיקה

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


השתמשתי בממשק ה- API של FCC מעט אבל עבור קואורדינטות של 19,000, התהליך ארך מספר שעות. נראה שהתסריט של GrantD71 עבד די מהר בכמה מאות, אך היו לי בעיות בהגדרת הסט שלי במועד. הצטרפות מרחבית של Geopandas היא הדרך המהירה ביותר שמצאתי עד כה. טפסי הצורה של TIGER נמצאים באתר המפקד. ייבא אותם ואת הקואורדינטות שלך כמו להלן (יתכן שתצטרך להקרין למערכת הפניה שונה של CRS).

ייבא גואפנדות כ- gpd ייבוא ​​פנדות כ- pd מ shapely.geometry ייבוא ​​נקודת ייבוא ​​os os.environ ['PROJ_LIB'] = r'C:  Users  root  Anaconda3  Library  share  proj 'census_tracts = gpd.read_file (r' C:  משתמשים  root  הורדות  tl_2018_06_bg  tl_2018_06_bg.shp ') points_df = pd.read_csv (r "C:  Users  root  Desktop  Desktop  airbnb_LAarea.csv", index_col = 0) גיאומטריה = [נקודה (xy ) עבור xy ב- zip (points_df.lng, points_df.lat)] crs = {'init': 'epsg: 4326'} gdf = gpd.GeoDataFrame (points_df, crs = crs, geometry = geometry) merged_file = gpd.sjoin ( gdf, census_tracts, how = "left", op = "within") merged_df = pd.DataFrame (merged_file)