البرمجة الجيومكانية

# إنشاء بيئة
conda create -n gis python=3.11 -y
conda activate gis

# حزم شائعة
conda install -c conda-forge geopandas pyogrio rtree shapely fiona -y
conda install -c conda-forge rasterio rioxarray xarray -y
pip install matplotlib contextily rich tqdm

ملاحظة: على Windows، تأكد من تثبيت Microsoft C++ Build Tools عند الحاجة لبعض الحزم.

import rasterio
from rasterio.plot import reshape_as_raster
import numpy as np

with rasterio.open("data/dem.tif") as src:
    dem = src.read(1)
    profile = src.profile
    print(profile["crs"], profile["transform"])
    print(np.nanmin(dem), np.nanmax(dem))

import rioxarray as rxr
import numpy as np

dem = rxr.open_rasterio("data/dem.tif").squeeze()
# تدرجات تقريبية (متر/بكسل)
dx = dem.differentiate("x")
dy = dem.differentiate("y")
slope = np.degrees(np.arctan(np.hypot(dx, dy)))

slope_da = dem.copy(data=slope)
slope_da.rio.to_raster("out/slope_deg.tif")

إذا أردت مؤشرات طبيعية مثل NDVI لاحقًا، ستجد نموذجًا جاهزًا في قسم المؤشرات أدناه ضمن بايثون.

ArcPy يعمل من بيئة ArcGIS Pro Python. مثال: قص طبقة، ثم Buffer، ثم حساب طول الطريق ضمن منطقة الدراسة.

import arcpy
arcpy.env.overwriteOutput = True

study = r"C:\gis\study.gdb\boundary"
roads = r"C:\gis\data.gdb\roads"
out_gdb = r"C:\gis\out.gdb"

clip_out = arcpy.analysis.Clip(roads, study, f"{out_gdb}\\roads_clip")
buf_out  = arcpy.analysis.Buffer(clip_out, f"{out_gdb}\\roads_buf_50m", "50 Meters")

arcpy.management.AddField(clip_out, "len_m", "DOUBLE")
arcpy.management.CalculateGeometryAttributes(clip_out,
    [["len_m","LENGTH_GEODESIC"]], length_unit="METERS")
print("done.")

نصيحة: احفظ التسلسل كـ ModelBuilder ثم صدّره ككود بايثون لتوثيق خطوات التحليل.

مقتطف موجز لاستخدام ArcGIS Maps SDK for JS (نسخة 4.x).

<link rel="stylesheet" href="https://js.arcgis.com/4.29/esri/themes/light/main.css">
<div id="view" style="height:400px"></div>
<script src="https://js.arcgis.com/4.29/"></script>
<script>
require(["esri/Map","esri/views/MapView","esri/layers/FeatureLayer"],(Map,MapView,FeatureLayer)=>{
  const map = new Map({basemap:"streets-navigation-vector"});
  const view = new MapView({container:"view", map, center:[46.7,24.7], zoom:6});
  const fl = new FeatureLayer({url:"YOUR/FeatureServer/0"});
  map.add(fl);
});
</script>

يمكن ضبط الرموز والـ popups والفلترة من خصائص FeatureLayer.

نموذج سكربت: يقرأ طبقة، يحسب طول الطرق لكل حي، يصدّر CSV وصورة بسيطة.

import geopandas as gpd
import matplotlib.pyplot as plt

districts = gpd.read_file("data/districts.gpkg").to_crs(epsg=32638)
roads     = gpd.read_file("data/roads.gpkg").to_crs(districts.crs)

length = roads.length.groupby(roads.sindex.query_bulk(districts.geometry, predicate="intersects")[1]).sum()
districts["roads_km"] = (length.reindex(range(len(districts))).fillna(0) / 1000).values

districts[["name","roads_km"]].to_csv("out/roads_by_district.csv", index=False)

ax = districts.plot(column="roads_km", cmap="viridis", legend=True, figsize=(6,6))
plt.tight_layout(); plt.savefig("out/roads_map.png", dpi=180)

شغّليه يدويًا، أو عبر جدولة النظام (Task Scheduler/cron) لتنفيذ يومي/أسبوعي.

CREATE DATABASE gisdb;
\c gisdb;
CREATE EXTENSION postgis;
-- تحقق
SELECT postgis_full_version();

ogr2ogr -f "PostgreSQL" PG:"dbname=gisdb user=postgres password=*" \
  roads.shp -nln public.roads -nlt PROMOTE_TO_MULTI -lco GEOMETRY_NAME=geom -lco FID=gid -overwrite

يفضل التخزين بـ SRID مناسب (UTM لمنطقتك) لسهولة القياسات بالأمتار.

1. تثبيت GeoServer (يفضل الإصدار المتوافق مع JRE حديث).
2. إضافة Store من نوع PostGIS وربطه بقاعدة gisdb.
3. إنشاء Layer لكل جدول مكاني وتعيين SRS صحيح.
4. اختبار WMS/WFS من تبويب Preview.

عند الاستهلاك من JavaScript، استخدمي عناوين WMS/WFS/WMTS الرسمية وقيود CORS حسب المضيف.

<sld:StyledLayerDescriptor version="1.0.0"
 xmlns:sld="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc">
  <sld:NamedLayer>
    <sld:Name>roads</sld:Name>
    <sld:UserStyle>
      <sld:FeatureTypeStyle>
        <sld:Rule>
          <sld:LineSymbolizer>
            <sld:Stroke>
              <sld:CssParameter name="stroke">#6b4dff</sld:CssParameter>
              <sld:CssParameter name="stroke-width">1.5</sld:CssParameter>
            </sld:Stroke>
          </sld:LineSymbolizer>
        </sld:Rule>
      </sld:FeatureTypeStyle>
    </sld:UserStyle>
  </sld:NamedLayer>
</sld:StyledLayerDescriptor>

import rasterio
from rasterio.windows import Window
from pathlib import Path
import numpy as np

src_path = "data/imagery.tif"
out_dir  = Path("out/tiles"); out_dir.mkdir(parents=True, exist_ok=True)

tile = 512
with rasterio.open(src_path) as src:
    for r in range(0, src.height, tile):
        for c in range(0, src.width, tile):
            w = Window(c, r, min(tile, src.width-c), min(tile, src.height-r))
            arr = src.read(window=w)
            if arr.size == 0: continue
            profile = src.profile
            profile.update({"height":arr.shape[1], "width":arr.shape[2],
                            "transform": rasterio.windows.transform(w, src.transform)})
            out = out_dir / f"tile_{r}_{c}.tif"
            with rasterio.open(out, "w", **profile) as dst:
                dst.write(arr)

بعدها استخدم أداة تعليم (Labeling) لرسم الأقنعة، ثم قسم إلى train/val.

# ملاحظة: هذا مخطط هيكلي مبسط (ليس جاهزًا للإنتاج)
import torch, torch.nn as nn, torch.optim as optim

class MiniUNet(nn.Module):
    def _init(self): super().init_()
    # طبقات مختصرة جداً للعرض فقط

# dataset/dataloader لقراءة المربعات (X: صور، y: أقنعة)
# loop: train for N epochs, save best model

# الحفظ:
# torch.save(model.state_dict(), "out/model.pt")

في المشاريع الفعلية يُفضّل استخدام مكتبات جاهزة لـ Remote Sensing مثل torchgeo أو segmentation-models-pytorch مع ضبط Augmentations ومقاييس التقييم.