AcadGIS — raster & vector layers (0.1.4)¶

A generic layer system: drop any raster (NDVI, land cover, climate, satellite) or any vector (your shapefile/GeoJSON, roads, cities) onto a map, plus a topography layer and XYZ basemaps. Every function is an add_*(ax, …) overlay, so it composes with plot() and every study_area panel.

function	what	needs
add_raster	GeoTIFF — continuous / categorical / RGB	acadgis[terrain]
add_layer	any vector, auto-styled, label control	core
add_basemap / add_satellite	XYZ tiles (6 styles)	acadgis[basemap]
add_topography	DEM hillshade + hypsometry	acadgis[terrain]
add_cities	Natural Earth places, ranked	core (+net)
add_roads	OpenStreetMap roads	core (+net)
add_landcover	ESA WorldCover 10 m (curated)	acadgis[terrain]+net
add_ndvi	Sentinel-2 NDVI (curated)	acadgis[terrain]+net

Z-order ladder: sea 0 · basemap/raster 1 · land 5 · rivers/roads 6 · points/labels 8. Label control (labels=): None hide · "one"/1 single · int top-N · "all" every feature.

In [1]:

import acadgis as agis
import tempfile, os
import numpy as np, rasterio                       # (demo-data generation only)
from rasterio.transform import from_bounds
print("AcadGIS", agis.__version__)

# --- base geometry ---
bd0 = agis.load_boundaries("Bangladesh", "country")
bd1 = agis.load_boundaries("Bangladesh", "division")
w, s, e, n = bd0.total_bounds

# --- DEMO DATA: synthesize an NDVI + a land-cover GeoTIFF over Bangladesh ---
# (In real use you skip this and pass your own .tif path to agis.add_raster.)
W = H = 480
LON, LAT = np.meshgrid(np.linspace(w, e, W), np.linspace(n, s, H))
ndvi = np.clip(0.55 + 0.25*np.cos(np.radians((LAT-23)*6)) + 0.20*np.sin(np.radians((LON-90)*5))
               + 0.18*np.exp(-(((LON-89.5)**2+(LAT-24.8)**2)/1.2))
               - 0.22*np.exp(-(((LON-90.4)**2+(LAT-23.8)**2)/0.25)), 0, 1).astype("float32")
tr = from_bounds(w, s, e, n, W, H)
NDVI_TIF = os.path.join(tempfile.gettempdir(), "demo_ndvi.tif")
with rasterio.open(NDVI_TIF, "w", driver="GTiff", height=H, width=W, count=1,
                   dtype="float32", crs="EPSG:4326", transform=tr, nodata=-9999) as d:
    d.write(ndvi, 1)

lc = np.full((H, W), 3, "uint8")                    # 3 cropland
lc[ndvi > 0.72] = 4                                 # 4 forest
lc[(ndvi > 0.45) & (ndvi <= 0.55)] = 5             # 5 wetland
lc[ndvi < 0.32] = 1                                 # 1 water
lc[np.exp(-(((LON-90.4)**2+(LAT-23.8)**2)/0.18)) > 0.5] = 2   # 2 built-up (Dhaka)
LC_TIF = os.path.join(tempfile.gettempdir(), "demo_landcover.tif")
with rasterio.open(LC_TIF, "w", driver="GTiff", height=H, width=W, count=1,
                   dtype="uint8", crs="EPSG:4326", transform=tr, nodata=0) as d:
    d.write(lc, 1)
LC_CLASSES = {1: ("#4a90d9", "Water"), 2: ("#d1495b", "Built-up"),
              3: ("#e9d8a6", "Cropland"), 4: ("#2a9d4a", "Forest"), 5: ("#76c893", "Wetland")}
print("demo rasters ready")

import acadgis as agis import tempfile, os import numpy as np, rasterio # (demo-data generation only) from rasterio.transform import from_bounds print("AcadGIS", agis.__version__) # --- base geometry --- bd0 = agis.load_boundaries("Bangladesh", "country") bd1 = agis.load_boundaries("Bangladesh", "division") w, s, e, n = bd0.total_bounds # --- DEMO DATA: synthesize an NDVI + a land-cover GeoTIFF over Bangladesh --- # (In real use you skip this and pass your own .tif path to agis.add_raster.) W = H = 480 LON, LAT = np.meshgrid(np.linspace(w, e, W), np.linspace(n, s, H)) ndvi = np.clip(0.55 + 0.25*np.cos(np.radians((LAT-23)*6)) + 0.20*np.sin(np.radians((LON-90)*5)) + 0.18*np.exp(-(((LON-89.5)**2+(LAT-24.8)**2)/1.2)) - 0.22*np.exp(-(((LON-90.4)**2+(LAT-23.8)**2)/0.25)), 0, 1).astype("float32") tr = from_bounds(w, s, e, n, W, H) NDVI_TIF = os.path.join(tempfile.gettempdir(), "demo_ndvi.tif") with rasterio.open(NDVI_TIF, "w", driver="GTiff", height=H, width=W, count=1, dtype="float32", crs="EPSG:4326", transform=tr, nodata=-9999) as d: d.write(ndvi, 1) lc = np.full((H, W), 3, "uint8") # 3 cropland lc[ndvi > 0.72] = 4 # 4 forest lc[(ndvi > 0.45) & (ndvi <= 0.55)] = 5 # 5 wetland lc[ndvi < 0.32] = 1 # 1 water lc[np.exp(-(((LON-90.4)**2+(LAT-23.8)**2)/0.18)) > 0.5] = 2 # 2 built-up (Dhaka) LC_TIF = os.path.join(tempfile.gettempdir(), "demo_landcover.tif") with rasterio.open(LC_TIF, "w", driver="GTiff", height=H, width=W, count=1, dtype="uint8", crs="EPSG:4326", transform=tr, nodata=0) as d: d.write(lc, 1) LC_CLASSES = {1: ("#4a90d9", "Water"), 2: ("#d1495b", "Built-up"), 3: ("#e9d8a6", "Cropland"), 4: ("#2a9d4a", "Forest"), 5: ("#76c893", "Wetland")} print("demo rasters ready")

AcadGIS 0.1.4
demo rasters ready

1 · add_raster — continuous (NDVI)¶

In [2]:

# Continuous raster — clip to the country, choose a colormap, add a colorbar.
fig, ax = agis.plt.subplots(figsize=(8, 9))
agis.add_raster(
    ax, NDVI_TIF,                      # any GeoTIFF / NetCDF path (or rioxarray DataArray)
    area=bd0,                          # clip to a geometry (GeoDataFrame/GeoSeries/shapely)
    cmap="YlGn", vmin=0, vmax=1,       # continuous colormap + range
    colorbar=True, colorbar_label="NDVI (synthetic)",
    opacity=1.0,                       # blend with layers beneath
)
agis.add_layer(ax, bd1, facecolor="none", edgecolor="#555", linewidth=0.6)   # outlines on top
ax.set_xlim(w, e); ax.set_ylim(s, n)
ax.set_title("add_raster — continuous (NDVI)", fontsize=12, fontweight="bold")
agis.show()

# Continuous raster — clip to the country, choose a colormap, add a colorbar. fig, ax = agis.plt.subplots(figsize=(8, 9)) agis.add_raster( ax, NDVI_TIF, # any GeoTIFF / NetCDF path (or rioxarray DataArray) area=bd0, # clip to a geometry (GeoDataFrame/GeoSeries/shapely) cmap="YlGn", vmin=0, vmax=1, # continuous colormap + range colorbar=True, colorbar_label="NDVI (synthetic)", opacity=1.0, # blend with layers beneath ) agis.add_layer(ax, bd1, facecolor="none", edgecolor="#555", linewidth=0.6) # outlines on top ax.set_xlim(w, e); ax.set_ylim(s, n) ax.set_title("add_raster — continuous (NDVI)", fontsize=12, fontweight="bold") agis.show()

2 · add_raster — categorical (land cover)¶

In [3]:

# Categorical raster — a class map + automatic legend (classes={value:(color,label)}).
fig, ax = agis.plt.subplots(figsize=(8, 9))
agis.add_raster(
    ax, LC_TIF, area=bd0,
    categorical=True, classes=LC_CLASSES,
    legend=True, legend_loc="lower left", legend_title="Land cover",
)
agis.add_layer(ax, bd1, facecolor="none", edgecolor="#333", linewidth=0.5)
ax.set_xlim(w, e); ax.set_ylim(s, n)
ax.set_title("add_raster — categorical (land cover)", fontsize=12, fontweight="bold")
agis.show()

# Categorical raster — a class map + automatic legend (classes={value:(color,label)}). fig, ax = agis.plt.subplots(figsize=(8, 9)) agis.add_raster( ax, LC_TIF, area=bd0, categorical=True, classes=LC_CLASSES, legend=True, legend_loc="lower left", legend_title="Land cover", ) agis.add_layer(ax, bd1, facecolor="none", edgecolor="#333", linewidth=0.5) ax.set_xlim(w, e); ax.set_ylim(s, n) ax.set_title("add_raster — categorical (land cover)", fontsize=12, fontweight="bold") agis.show()

3 · add_basemap (satellite) + add_cities¶

In [4]:

# Real satellite basemap + cities, with label control.
#   labels: None/False = hide · True/"all" = every city · "one"/1 = single · int N = top-N
dist = agis.load_boundaries("Bangladesh", "district", within="Dhaka")
gz = dist[dist["NAME_2"] == "Gazipur"]
gw, gs, ge, gn = gz.total_bounds
fig, ax = agis.plt.subplots(figsize=(8.5, 8))
ax.set_xlim(gw - .05, ge + .05); ax.set_ylim(gs - .05, gn + .05)
agis.add_basemap(ax, style="satellite")            # satellite · osm · light · dark · terrain · toner
agis.add_layer(ax, gz, facecolor="none", edgecolor="#ffd000", linewidth=2.0)
agis.add_cities(ax, top=12, labels="all", color="#fff70a",
                label_color="white", markersize=30)   # NE populated places, ranked by pop
ax.set_title("add_basemap (satellite) + add_cities", fontsize=12, fontweight="bold")
agis.show()

# Real satellite basemap + cities, with label control. # labels: None/False = hide · True/"all" = every city · "one"/1 = single · int N = top-N dist = agis.load_boundaries("Bangladesh", "district", within="Dhaka") gz = dist[dist["NAME_2"] == "Gazipur"] gw, gs, ge, gn = gz.total_bounds fig, ax = agis.plt.subplots(figsize=(8.5, 8)) ax.set_xlim(gw - .05, ge + .05); ax.set_ylim(gs - .05, gn + .05) agis.add_basemap(ax, style="satellite") # satellite · osm · light · dark · terrain · toner agis.add_layer(ax, gz, facecolor="none", edgecolor="#ffd000", linewidth=2.0) agis.add_cities(ax, top=12, labels="all", color="#fff70a", label_color="white", markersize=30) # NE populated places, ranked by pop ax.set_title("add_basemap (satellite) + add_cities", fontsize=12, fontweight="bold") agis.show()

4 · add_layer — any vector (lines + labelled points)¶

In [5]:

# Generic vector — bring any GeoDataFrame/file. Auto-styles polygon/line/point.
rivers = agis.load_rivers(scale="10m")
sites = agis.gpd.GeoDataFrame(
    {"name": [f"Site {i+1}" for i in range(5)],
     "geometry": [agis.gpd.points_from_xy([90.0, 91.2, 89.3, 90.8, 88.9],
                                           [24.5, 23.2, 22.9, 25.3, 24.0])[i] for i in range(5)]},
    crs="EPSG:4326")
fig, ax = agis.plt.subplots(figsize=(8.5, 9))
agis.plot(bd1, ax=ax, palette="pastel", graticule=False, north_arrow=False, scale_bar=False)
agis.add_layer(ax, rivers, area=bd0, color="#2b7bba", linewidth=0.9, label="Rivers")  # lines
agis.add_layer(ax, sites, color="#c0392b", markersize=55, marker="^",
               labels="all", label_field="name", label_size=8,            # point labels
               label="Field sites", legend=True)
ax.set_xlim(w, e); ax.set_ylim(s, n)
ax.set_title("add_layer — rivers (lines) + sites (points, labelled)", fontsize=12, fontweight="bold")
agis.show()

# Generic vector — bring any GeoDataFrame/file. Auto-styles polygon/line/point. rivers = agis.load_rivers(scale="10m") sites = agis.gpd.GeoDataFrame( {"name": [f"Site {i+1}" for i in range(5)], "geometry": [agis.gpd.points_from_xy([90.0, 91.2, 89.3, 90.8, 88.9], [24.5, 23.2, 22.9, 25.3, 24.0])[i] for i in range(5)]}, crs="EPSG:4326") fig, ax = agis.plt.subplots(figsize=(8.5, 9)) agis.plot(bd1, ax=ax, palette="pastel", graticule=False, north_arrow=False, scale_bar=False) agis.add_layer(ax, rivers, area=bd0, color="#2b7bba", linewidth=0.9, label="Rivers") # lines agis.add_layer(ax, sites, color="#c0392b", markersize=55, marker="^", labels="all", label_field="name", label_size=8, # point labels label="Field sites", legend=True) ax.set_xlim(w, e); ax.set_ylim(s, n) ax.set_title("add_layer — rivers (lines) + sites (points, labelled)", fontsize=12, fontweight="bold") agis.show()

5 · add_topography — shaded relief from a DEM¶

In [6]:

# Topography — downloads a DEM and renders hypsometric + hillshaded relief.
fig, ax = agis.plt.subplots(figsize=(8, 8.5))
agis.add_topography(
    ax, area=(91.7, 21.9, 92.5, 23.2),     # Chittagong Hill Tracts (name / bbox / GeoDataFrame)
    cmap="terrain", hillshade=True,
    colorbar=True, colorbar_label="Elevation (m)",
    max_size=900,                          # DEM resolution cap (speed vs detail)
)
ax.set_title("add_topography — shaded relief (DEM)", fontsize=12, fontweight="bold")
agis.show()

# Topography — downloads a DEM and renders hypsometric + hillshaded relief. fig, ax = agis.plt.subplots(figsize=(8, 8.5)) agis.add_topography( ax, area=(91.7, 21.9, 92.5, 23.2), # Chittagong Hill Tracts (name / bbox / GeoDataFrame) cmap="terrain", hillshade=True, colorbar=True, colorbar_label="Elevation (m)", max_size=900, # DEM resolution cap (speed vs detail) ) ax.set_title("add_topography — shaded relief (DEM)", fontsize=12, fontweight="bold") agis.show()

6 · add_roads — OpenStreetMap roads over a basemap¶

In [7]:

# Roads — OpenStreetMap road network over a satellite basemap (widths by class).
fig, ax = agis.plt.subplots(figsize=(8.5, 8))
ax.set_xlim(90.33, 90.46); ax.set_ylim(23.72, 23.85)        # central Dhaka
agis.add_basemap(ax, style="light")
agis.add_roads(ax, classes=("motorway", "trunk", "primary", "secondary"),
               color="#d1495b", by_class=True)               # auto width per highway class
agis.add_cities(ax, top=1, labels="one", color="#111", label_color="#111", markersize=40)
ax.set_title("add_roads (OSM) over a basemap — Dhaka", fontsize=12, fontweight="bold")
agis.show()

# Roads — OpenStreetMap road network over a satellite basemap (widths by class). fig, ax = agis.plt.subplots(figsize=(8.5, 8)) ax.set_xlim(90.33, 90.46); ax.set_ylim(23.72, 23.85) # central Dhaka agis.add_basemap(ax, style="light") agis.add_roads(ax, classes=("motorway", "trunk", "primary", "secondary"), color="#d1495b", by_class=True) # auto width per highway class agis.add_cities(ax, top=1, labels="one", color="#111", label_color="#111", markersize=40) ax.set_title("add_roads (OSM) over a basemap — Dhaka", fontsize=12, fontweight="bold") agis.show()

7 · Compose a full figure¶

In [8]:

# Compose a research-grade figure — basemap + boundaries + rivers + ranked cities.
fig, ax = agis.plt.subplots(figsize=(9, 10))
ax.set_xlim(w, e); ax.set_ylim(s, n)
agis.add_basemap(ax, style="light", alpha=0.9)             # z 0
agis.add_layer(ax, bd1, facecolor="none", edgecolor="#444", linewidth=0.8)   # z 5
agis.add_layer(ax, agis.load_rivers(scale="10m"), area=bd0, color="#2b7bba", linewidth=0.7)  # z 6
agis.add_cities(ax, top=8, labels=5, color="#111", label_color="#111",       # show 8, label top 5
                markersize=34, marker="o")                  # z 8
ax.set_title("Composed layers — basemap · boundaries · rivers · cities",
             fontsize=12, fontweight="bold")
agis.show()

# Compose a research-grade figure — basemap + boundaries + rivers + ranked cities. fig, ax = agis.plt.subplots(figsize=(9, 10)) ax.set_xlim(w, e); ax.set_ylim(s, n) agis.add_basemap(ax, style="light", alpha=0.9) # z 0 agis.add_layer(ax, bd1, facecolor="none", edgecolor="#444", linewidth=0.8) # z 5 agis.add_layer(ax, agis.load_rivers(scale="10m"), area=bd0, color="#2b7bba", linewidth=0.7) # z 6 agis.add_cities(ax, top=8, labels=5, color="#111", label_color="#111", # show 8, label top 5 markersize=34, marker="o") # z 8 ax.set_title("Composed layers — basemap · boundaries · rivers · cities", fontsize=12, fontweight="bold") agis.show()

8 · add_landcover — curated ESA WorldCover (real data)¶

In [9]:

# Curated land cover — ESA WorldCover 10 m fetched automatically for the area.
AREA = (90.28, 23.62, 90.55, 23.95)        # bbox (minx,miny,maxx,maxy) — or a name / GeoDataFrame
dhk = agis.load_boundaries("Bangladesh", "district", within="Dhaka")
fig, ax = agis.plt.subplots(figsize=(8.5, 8))
ax.set_xlim(AREA[0], AREA[2]); ax.set_ylim(AREA[1], AREA[3])
agis.add_landcover(ax, AREA, year=2021)    # real 10 m land cover + official 11-class legend
agis.add_layer(ax, dhk, facecolor="none", edgecolor="#111", linewidth=0.7)
ax.set_title("add_landcover — ESA WorldCover 10 m (real data)", fontsize=12, fontweight="bold")
agis.show()

# Curated land cover — ESA WorldCover 10 m fetched automatically for the area. AREA = (90.28, 23.62, 90.55, 23.95) # bbox (minx,miny,maxx,maxy) — or a name / GeoDataFrame dhk = agis.load_boundaries("Bangladesh", "district", within="Dhaka") fig, ax = agis.plt.subplots(figsize=(8.5, 8)) ax.set_xlim(AREA[0], AREA[2]); ax.set_ylim(AREA[1], AREA[3]) agis.add_landcover(ax, AREA, year=2021) # real 10 m land cover + official 11-class legend agis.add_layer(ax, dhk, facecolor="none", edgecolor="#111", linewidth=0.7) ax.set_title("add_landcover — ESA WorldCover 10 m (real data)", fontsize=12, fontweight="bold") agis.show()

9 · add_ndvi — curated Sentinel-2 NDVI (real data)¶

In [10]:

# Curated NDVI — least-cloudy Sentinel-2 scene in the window, computed live.
fig, ax = agis.plt.subplots(figsize=(8.5, 8))
ax.set_xlim(AREA[0], AREA[2]); ax.set_ylim(AREA[1], AREA[3])
agis.add_ndvi(ax, AREA, start="2024-01-01", end="2024-04-01",   # widen for cloudy regions
              max_cloud=8, cmap="RdYlGn")
agis.add_layer(ax, dhk, facecolor="none", edgecolor="#111", linewidth=0.7)
ax.set_title(f"add_ndvi — Sentinel-2 NDVI ({ax._acadgis_ndvi_date}, real)",
             fontsize=12, fontweight="bold")
agis.show()

# Curated NDVI — least-cloudy Sentinel-2 scene in the window, computed live. fig, ax = agis.plt.subplots(figsize=(8.5, 8)) ax.set_xlim(AREA[0], AREA[2]); ax.set_ylim(AREA[1], AREA[3]) agis.add_ndvi(ax, AREA, start="2024-01-01", end="2024-04-01", # widen for cloudy regions max_cloud=8, cmap="RdYlGn") agis.add_layer(ax, dhk, facecolor="none", edgecolor="#111", linewidth=0.7) ax.set_title(f"add_ndvi — Sentinel-2 NDVI ({ax._acadgis_ndvi_date}, real)", fontsize=12, fontweight="bold") agis.show()