تحليل البيانات المكانية في R: من الحدود والخرائط إلى موران والنقاط الساخنة

لديك جدول مؤشرات على مستوى البلديات أو المحافظات: معدل تغطية خدمة، أو عدد المراكز الصحية لكل ألف نسمة، أو نسبة الأسر تحت خط معين. الإدارة لا تسأل عن المتوسط العام، بل تسأل: أين تتركز الفجوة؟ وهل المناطق المحرومة متجاورة فتحتاج تدخلًا إقليميًا، أم متناثرة فتحتاج حلولًا محلية؟ خريطة ملونة في Excel لا تجيب عن هذا السؤال، لأنها لا تختبر إن كان النمط الذي تراه أقوى مما تنتجه الصدفة. هنا يأتي دور تحليل البيانات المكانية في R: قراءة الحدود الإدارية بإحداثياتها الصحيحة، وربطها بالجداول، ورسم خرائط موضوعية لا تضلل، ثم قياس الارتباط الذاتي المكاني وتحديد النقاط الساخنة باختبار إحصائي لا بالتخمين. نستخدم في الأمثلة بيانات متلازمة موت الرضع المفاجئ (SIDS) لمقاطعات ولاية كارولاينا الشمالية التي تأتي مع حزمة sf، فيمكنك تشغيل كل سطر على جهازك، والنتائج التي ستظهر لك توضيحية لهذه البيانات فقط.

البيانات المتجهية والبيانات الشبكية

أول قرار في أي مشروع مكاني هو نوع البيانات. البيانات المتجهية (Vector) تمثل العالم بأشكال هندسية: نقاط للمدارس والمراكز الصحية، وخطوط للطرق وشبكات المياه، ومضلعات للحدود الإدارية والأحياء. أما البيانات الشبكية (Raster) فتقسم المساحة إلى خلايا منتظمة، لكل خلية قيمة: ارتفاع، أو حرارة، أو غطاء أرضي، أو كثافة سكان تقديرية. المؤشرات الإدارية والإحصائية تعيش عادةً في العالم المتجهي، بينما تأتي صور الأقمار الصناعية والنماذج البيئية شبكية.

كثير من المشاريع تجمع النوعين: تحسب مثلًا متوسط كثافة السكان الشبكية داخل كل حي باستخراج القيم بالمضلعات (الدالة extract() في terra)، ثم تكمل التحليل على المستوى المتجهي. يركز هذا المقال على البيانات المتجهية لأنها لغة المؤشرات والتخطيط.

أنظمة الإحداثيات المرجعية: مصدر أكثر الأخطاء صمتًا

كل طبقة مكانية معرفة بنظام إحداثيات مرجعي (CRS). الأنظمة الجغرافية، مثل WGS 84 (رمزه EPSG:4326)، تعبر عن المواقع بخطوط الطول والعرض بالدرجات، والدرجة ليست وحدة طول ثابتة: طول درجة خط الطول يتقلص كلما ابتعدت عن خط الاستواء. أما الأنظمة المسقطة (Projected)، مثل نطاق UTM المناسب لمنطقتك، فتحول سطح الأرض إلى مستوى بوحدات مترية، فتصبح المساحات والمسافات والحواجز (Buffers) ذات معنى مباشر.

قاعدتان عمليتان: لا تجمع طبقتين بنظامين مختلفين قبل التحويل، واحسب المساحات والمسافات في نظام مسقط مناسب للمنطقة. تستخدم sf في إصداراتها الحديثة مكتبة s2 للحساب على سطح الكرة حين تكون الإحداثيات جغرافية، فتعطي مساحات معقولة حتى بالدرجات، لكن العمل في نظام مسقط يبقى أوضح وأسرع في عمليات الحواجز والتقاطع.

library(sf)

nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
st_crs(nc)$epsg            # 4267: NAD27, geographic (degrees)
st_geometry_type(nc)[1]    # MULTIPOLYGON

# Project to NAD83 / North Carolina state plane (metres) before measuring
nc_m <- st_transform(nc, 32119)
nc_m$area_km2 <- as.numeric(st_area(nc_m)) / 1e6
head(st_drop_geometry(nc_m)[, c("NAME", "FIPS", "BIR79", "SID79", "area_km2")])

البيانات مخزنة بنظام NAD27 الجغرافي، فحولناها إلى نظام الولاية المسقط بالأمتار قبل حساب المساحة. في بيانات المنطقة العربية ستجد غالبًا ملفات بنظام WGS 84، أو بأنظمة وطنية محلية، أو ملفات Shapefile بلا ملف .prj أصلًا. في الحالة الأخيرة لا تخمن: اسأل الجهة المصدرة عن النظام، ثم عينه بالدالة st_set_crs()، وبعدها فقط حول بالدالة st_transform(). الأولى تكتب اسم النظام على البيانات دون تغيير الإحداثيات، والثانية تحسب إحداثيات جديدة، والخلط بينهما يزيح الطبقة كلها.

قراءة الحدود وربطها بالجداول

ستحتاج باستمرار إلى نوعين من الربط. الربط الوصفي (Attribute Join) يضيف جدول المؤشرات إلى طبقة الحدود عبر مفتاح مشترك. والربط المكاني (Spatial Join) يستخدم الهندسة نفسها: أي نقطة تقع في أي مضلع، وأي حي يتقاطع مع أي منطقة خدمة.

المفتاح في الربط الوصفي يجب أن يكون رمزًا إداريًا ثابتًا، لا اسمًا. أسماء المحافظات والبلديات بالعربية تكتب بصيغ متعددة: بالهمزة أو دونها، بالتاء المربوطة أو الهاء، مع «ال» أو دونها. وكل اختلاف إملائي يتحول إلى مضلع فارغ على الخريطة دون أي رسالة خطأ. إذا لم يتوفر رمز موحد، فابن جدول مطابقة مرة واحدة واحفظه مع المشروع ليستخدم في كل تحديث.

library(dplyr)

# Attribute join: an indicator table keyed by the county FIPS code (synthetic values)
set.seed(1)
indicators <- data.frame(FIPS = nc$FIPS, clinics = rpois(nrow(nc), 3))
nc_j <- left_join(nc_m, indicators, by = "FIPS")
stopifnot(!anyNA(nc_j$clinics))     # every polygon must find its row

# Spatial join: synthetic facility points, counted per county
pts <- st_sf(id = 1:300, geometry = st_sample(nc_m, 300))
counts <- st_join(pts, nc_m["FIPS"]) |>
  st_drop_geometry() |>
  count(FIPS, name = "n_facilities")

nc_j <- left_join(nc_j, counts, by = "FIPS") |>
  mutate(n_facilities = coalesce(n_facilities, 0L))

سطر stopifnot ليس ترفًا: هو ما يكشف أن ثلاث بلديات لم تطابق قبل أن تظهر بيضاء في تقرير رسمي. ولاحظ أن الربط المكاني يتطلب أن تكون الطبقتان بالنظام المرجعي نفسه، وإلا أوقفتك sf برسالة خطأ صريحة، وهذا من حسن حظك.

الخرائط الموضوعية: معدلات لا أعداد

الخريطة الموضوعية (Choropleth) تلون كل وحدة إدارية حسب قيمة المؤشر. خطؤها الأشهر رسم الأعداد المطلقة: الوحدات الأكبر سكانًا ستظهر دائمًا «أسوأ» لأنها أكبر فقط. ارسم المعدلات أو النسب، أي العدد مقسومًا على السكان المعرضين للحدث. والقرار الثاني هو طريقة التصنيف: الفئات المتساوية المدى تبرز القيم المتطرفة، والمئينات (Quantiles) توزع الوحدات بالتساوي على الألوان، والفواصل الطبيعية (Jenks) تبحث عن فجوات في التوزيع. لا توجد طريقة صحيحة دائمًا، لكن اذكر الطريقة في مفتاح الخريطة، لأن الخريطة نفسها قد تبدو مختلفة كليًا بتغيير التصنيف.

library(ggplot2)

nc_m$sids_rate <- 1000 * nc_m$SID79 / nc_m$BIR79      # deaths per 1,000 live births
brks <- unique(quantile(nc_m$sids_rate, probs = seq(0, 1, 0.2)))
nc_m$rate_q <- cut(nc_m$sids_rate, brks, include.lowest = TRUE)

ggplot(nc_m) +
  geom_sf(aes(fill = rate_q), colour = "white", linewidth = 0.1) +
  scale_fill_brewer(palette = "YlOrRd", name = "SIDS per 1,000n(quintiles)") +
  labs(title = "SIDS rate by county, 1979-84") +
  theme_void()

إن كنت تفضل tmap، فالسطر tm_shape(nc_m) + tm_polygons("sids_rate") يعطيك خريطة سريعة، وtmap_mode("view") يحولها إلى خريطة تفاعلية قابلة للتكبير. انتبه إلى أن الإصدار الرابع من tmap غير أسماء كثير من الوسائط (مثل fill وfill.scale)، فالأمثلة المكتوبة للإصدار الثالث على الإنترنت قد تعطي تحذيرات أو نتائج مختلفة.

وهناك مشكلة تخص المعدلات نفسها: الوحدات قليلة السكان تعطي معدلات متطرفة بسبب صغر المقام لا بسبب واقع مختلف، فحالة واحدة في بلدة صغيرة قد تضعها في أعلى فئة. الحل الشائع هو التمهيد البيزي التجريبي (Empirical Bayes) الذي يسحب المعدلات غير المستقرة نحو المتوسط بقدر عدم استقرارها، وتوفره حزمة spdep بالدالة EBest().

الارتباط الذاتي المكاني: هل النمط أقوى من الصدفة؟

يقوم التحليل المكاني على فكرة أن الأشياء المتقاربة أكثر تشابهًا من المتباعدة، وهي ما يعرف بقانون توبلر الأول في الجغرافيا. الارتباط الذاتي المكاني (Spatial Autocorrelation) يقيس هذا التشابه رقميًا. وقبل أي مقياس عليك أن تعرف «الجوار»، أي تحديد الوحدات التي تعد جارة لكل وحدة، ثم تحول الجوار إلى مصفوفة أوزان مكانية.

موران العام (Global Moran’s I)

يلخص مؤشر موران العام النمط في رقم واحد. القيمة الموجبة تعني أن القيم المرتفعة تجاور المرتفعة والمنخفضة تجاور المنخفضة، أي تكتلًا. والقيمة السالبة تعني تناوبًا يشبه رقعة الشطرنج. والقيمة القريبة من توقعها تحت فرضية العشوائية، وهو -1/(n-1)، تعني غياب نمط مكاني واضح.

library(spdep)

nb <- poly2nb(nc_m, queen = TRUE)                     # contiguity neighbours
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)  # row-standardised weights

moran.test(nc_m$sids_rate, lw, zero.policy = TRUE)    # analytical test

set.seed(42)
moran.mc(nc_m$sids_rate, lw, nsim = 999, zero.policy = TRUE)  # permutation test

يعطيك moran.test قيمة I مع اختبار يعتمد على تقريب نظري، بينما يوزع moran.mc القيم عشوائيًا على الخريطة 999 مرة ويقارن القيمة الفعلية بهذه التوزيعات، وهو أكثر متانة حين يكون المؤشر ملتويًا كما هو حال المعدلات غالبًا. والأوزان المعيارية بالصفوف (style = "W") تجعل «القيمة المتأخرة مكانيًا» لكل وحدة متوسط قيم جيرانها، وهو التفسير الأسهل في التقارير.

موران المحلي (LISA) وخريطة التكتلات

موران العام يقول إن هناك تكتلًا، ولا يقول أين. المؤشرات المحلية للارتباط المكاني (LISA) التي اقترحها أنسلين تحسب مؤشرًا لكل وحدة، وتصنف الوحدات الدالة إلى أربع فئات: مرتفع محاط بمرتفع (High-High)، ومنخفض محاط بمنخفض (Low-Low)، وحالتان شاذتان مكانيًا: مرتفع وسط منخفض، ومنخفض وسط مرتفع.

lisa <- localmoran(nc_m$sids_rate, lw, zero.policy = TRUE)
p    <- lisa[, 5]                                   # two-sided p-value (5th column)
quad <- as.character(attr(lisa, "quadr")$mean)      # High-High, Low-Low, ...

nc_m$lisa <- ifelse(p < 0.05, quad, "Not significant")
table(nc_m$lisa)

# Stricter version: control the false discovery rate across 100 local tests
nc_m$lisa_fdr <- ifelse(p.adjust(p, method = "fdr") < 0.05, quad, "Not significant")

تنبيهان مهمان. الأول أنك تجري اختبارًا لكل وحدة، فمع مئة مقاطعة ستظهر بعض النتائج الدالة بالصدفة وحدها، ولهذا أضفنا نسخة مصححة بمعدل الاكتشاف الخاطئ؛ اعرض الاثنتين أو اشرح اختيارك. الثاني أن الاختبارات المحلية متداخلة لأن الجيران مشتركون بين الوحدات، لذلك يفضل كثيرون الاستدلال بالتبديل الشرطي، وتوفره الدالة localmoran_perm() في الإصدارات الحديثة من spdep.

النقاط الساخنة بإحصاءة Getis-Ord

إذا كان سؤال الإدارة مباشرًا، «أين تتركز القيم المرتفعة؟»، فإحصاءة Gi* لغيتس وأورد هي الأداة المعتادة. هي تقارن مجموع القيم في الوحدة وجيرانها بالمجموع المتوقع، وتعطي درجة معيارية z: الموجبة الكبيرة نقطة ساخنة، والسالبة الكبيرة نقطة باردة. النسخة النجمية تضم الوحدة نفسها إلى جوارها، ولهذا نستخدم include.self().

lw_self <- nb2listw(include.self(nb), style = "B")   # Gi* includes the unit itself
nc_m$gi_z <- as.numeric(localG(nc_m$sids_rate, lw_self))

nc_m$hotspot <- cut(nc_m$gi_z, breaks = c(-Inf, -1.96, 1.96, Inf),
                    labels = c("Cold spot", "Not significant", "Hot spot"))
table(nc_m$hotspot)

الحد 1.96 يقابل مستوى دلالة 0.05 لاختبار منفرد ثنائي الطرف، ولا يصحح للمقارنات المتعددة، فاعتبره عتبة استكشافية. وتحت هذه النتائج جميعًا افتراض أن مصفوفة الأوزان معقولة؛ جرب تعريفين للجوار على الأقل، فإذا تغيرت الخريطة جذريًا فهذا بحد ذاته نتيجة تستحق الذكر في التقرير.

الانحدار المكاني بإيجاز

حين تبني نموذج انحدار عاديًا على وحدات مكانية، كثيرًا ما تبقى البواقي مترابطة مكانيًا، فتصبح الأخطاء المعيارية متفائلة والاستنتاجات أقوى مما تسمح به البيانات. الخطوة الأولى اختبار موران على البواقي. فإن كان دالًا فأمامك نموذجان أساسيان: نموذج التأخر المكاني (Spatial Lag) الذي يفترض أن قيمة المتغير التابع في منطقة تتأثر بقيمه في المناطق المجاورة، كأثر انتشار أو تنافس على الخدمة، ونموذج الخطأ المكاني (Spatial Error) الذي يفترض أن الترابط في البواقي ناتج عن متغيرات غير مقاسة لها نمط مكاني. ولمراجعة أساسيات النموذج الخطي وافتراضاته اقرأ مقال تحليل الانحدار.

library(spatialreg)

nc_m$nwbir_pct <- 100 * nc_m$NWBIR79 / nc_m$BIR79
ols <- lm(sids_rate ~ nwbir_pct, data = nc_m)
lm.morantest(ols, lw, zero.policy = TRUE)      # spatial autocorrelation in residuals?

sar <- lagsarlm(sids_rate ~ nwbir_pct, data = nc_m, listw = lw, zero.policy = TRUE)
sem <- errorsarlm(sids_rate ~ nwbir_pct, data = nc_m, listw = lw, zero.policy = TRUE)
summary(sar)                                   # rho = spatial lag parameter
AIC(ols, sar, sem)
impacts(sar, listw = lw)                       # direct, indirect and total effects

في نموذج التأخر المكاني لا تفسر المعاملات مباشرةً كأثر حدي كما في الانحدار العادي، لأن التغير في منطقة ينتقل إلى جيرانها ثم يرتد إليها، ولهذا نستخدم impacts() للحصول على الآثار المباشرة وغير المباشرة والكلية. وللاختيار بين النموذجين تستخدم عادةً اختبارات مضاعف لاغرانج على بواقي النموذج العادي إلى جانب المعرفة النظرية بالظاهرة، فلا تختر النموذج بمعيار AIC وحده.

تطبيقات في التخطيط وتوزيع الخدمات

تتحول هذه الأدوات إلى قرارات حين تربطها بسؤال تخطيطي محدد. من الاستخدامات المعتادة في الوزارات والبلديات والمنظمات:

• إمكانية الوصول: حساب المسافة من كل تجمع سكاني إلى أقرب مدرسة أو مركز صحي، وتحديد السكان خارج نطاق مقبول تحدده الجهة المختصة.
• مناطق الخدمة: حواجز حول المرافق بنصف قطر يمثل نطاق الخدمة، ثم حساب السكان المغطين وغير المغطين.
• أولوية التدخل: مطابقة خريطة النقاط الساخنة للاحتياج مع خريطة توزيع الخدمات، لتحديد المناطق التي يجتمع فيها احتياج مرتفع وتغطية منخفضة.
• اختيار مواقع جديدة: مقارنة المواقع المرشحة بعدد السكان الذين يضيفهم كل منها إلى نطاق الخدمة.
• ربط المسوح بالجغرافيا: نقل تقديرات المسوح إلى الخرائط، مع الحذر من أن التقديرات المباشرة على المستوى المحلي قد تكون غير مستقرة، وهو ما يعالجه تقدير المناطق الصغيرة الذي نشرحه في مقال تحليل المسوح المعقدة وأوزانها.

مثال مختصر لحساب المسافة من مركز كل مقاطعة إلى أقرب مرفق من النقاط الاصطناعية التي أنشأناها سابقًا:

cent    <- st_centroid(st_geometry(nc_m))
nearest <- st_nearest_feature(cent, pts)          # index of the closest facility
nc_m$dist_km <- as.numeric(
  st_distance(cent, pts[nearest, ], by_element = TRUE)) / 1000
summary(nc_m$dist_km)

مركز المضلع الهندسي تقريب خشن: في الوحدات الكبيرة قد يقع في منطقة خالية من السكان. الأدق استخدام مراكز موزونة بالسكان، أو حساب المسافة لخلايا سكانية شبكية ثم تجميعها، وهنا يلتقي العالمان المتجهي والشبكي مرة أخرى. والمسافة الخطية نفسها تقريب، فحين يهم زمن الوصول الفعلي تحتاج إلى شبكة الطرق وتحليل الشبكات.

أخطاء شائعة

• دمج طبقات بأنظمة إحداثيات مختلفة، أو حساب المساحات والمسافات بالدرجات.
• الخلط بين st_set_crs() وst_transform().
• الربط بالأسماء العربية بدل الرموز الإدارية، دون فحص الوحدات التي لم تطابق.
• رسم أعداد مطلقة في خريطة موضوعية بدل المعدلات، أو عرض معدلات وحدات صغيرة جدًا دون تمهيد.
• تجاهل مشكلة الوحدة المساحية القابلة للتعديل (MAUP): النتائج قد تتغير بتغير مستوى التجميع أو رسم حدود الوحدات.
• تعميم نتائج المناطق على الأفراد، وهي المغالطة البيئية (Ecological Fallacy).
• قراءة كل وحدة دالة في LISA نتيجة مؤكدة دون تصحيح للمقارنات المتعددة.
• استخدام zero.policy = TRUE دون معرفة عدد الوحدات التي لا جيران لها ولماذا.
• الاعتماد على Shapefile في المشاريع الجديدة رغم قيوده، كأسماء الأعمدة القصيرة ومشكلات ترميز النصوص العربية، بينما GeoPackage ملف واحد أنظف.

أسئلة شائعة

هل أحتاج إلى R إذا كنت أستخدم QGIS أو ArcGIS؟

برامج نظم المعلومات الجغرافية ممتازة للتحرير والرسم والاستكشاف. ميزة R أن كل خطوة مكتوبة في نص برمجي يعاد تشغيله حين تصل بيانات جديدة، وأن الإحصاء المكاني والنمذجة متاحان في البيئة نفسها. كثير من الفرق تستخدم الاثنين: QGIS للمعاينة والإخراج، وR للتحليل القابل للتكرار.

أي مصفوفة أوزان أختار؟

للحدود الإدارية المتلاصقة يبدأ أغلب المحللين بالتجاور الملكي مع أوزان معيارية بالصفوف. وإذا كانت الوحدات متفاوتة الحجم جدًا أو فيها جزر، فأقرب k جيران أنسب. الأهم أن تبرر اختيارك نظريًا، وأن تختبر حساسية النتائج لتعريف بديل.

هل يعني موران الدال أن هناك سببًا مكانيًا؟

لا. هو يثبت أن القيم متكتلة أكثر مما تنتجه الصدفة، لكن التكتل قد ينتج عن متغيرات أخرى متكتلة بدورها، كالدخل أو الكثافة السكانية. التفسير السببي يحتاج إلى نموذج وتصميم مناسبين، والانحدار المكاني خطوة في هذا الاتجاه لا برهان نهائي.

هل يمكن تنفيذ التحليل نفسه في Python؟

نعم. تقابل geopandas حزمة sf في القراءة والربط والتحويل، وتوفر مكتبات PySAL، مثل libpysal وesda، مصفوفات الأوزان ومؤشرات موران العامة والمحلية. الاختيار بين اللغتين يعتمد على فريقك وبقية سلسلة العمل، وقد فصلنا ذلك في مقال Python أم R لتحليل البيانات.

التحليل المكاني الجيد يبدأ بنظام إحداثيات صحيح وربط مفحوص، ويمر بخرائط معدلات تذكر طريقة تصنيفها، وينتهي باختبارات تفصل النمط الحقيقي عن الصدفة. وإذا كانت لدى مؤسستك مؤشرات على مستوى المناطق وتحتاج إلى خرائط وتحليل تكتلات يدعم قرارات توزيع الخدمات، فيمكن لفريق مركز المستشار تنفيذ ذلك بلغتي R وPython ضمن خدمة تحليل البيانات المتقدم، والاستشارة الأولى مجانية.