Where Morocco’s jobs are missing

Why this map

The 2024 General Census of Population and Housing (RGPH 2024), run by the Haut-Commissariat au Plan (HCP) in September 2024, put Morocco’s legal population at 36,828,330 and its census unemployment rate at 21.3 percent (HCP, 2024). The national figure hides a sharp regional split. Casablanca-Settat holds about one Moroccan in five and sits below the national rate. The Oriental and Guelmim-Oued Noun, two peripheral regions, sit above 30 percent.

This notebook maps that split across the 12 regions created by the 2015 regionalization reform. It pulls the boundaries from an open API, joins them to HCP census figures, normalizes population by land area and checks each regional capital against the polygon it falls in.

Note

The census unemployment rate (21.3 percent) is higher than the rate from the HCP labour force survey (13.3 percent for 2024). The two instruments define job search differently, and the census counts every household. The census rate is the only one HCP publishes for all 12 regions individually; the labour force survey groups the three southern regions into one unit (HCP, 2025).

Data sources

The pipeline has three inputs.

• The regional polygons come from the geoBoundaries API (Runfola et al., 2020). The API is open and needs no key. A request returns a JSON record with the download link for the Morocco ADM1 file, and the notebook saves that file to data/ so that later runs work offline.
• Population and unemployment by region come from HCP’s RGPH 2024 releases and sit in data/regions_rgph2024.csv, keyed by ISO 3166-2 code. HCP has no public API, so these 12 rows were typed from its published tables.
• The 12 regional capitals are in data/regional_capitals.csv as latitude and longitude, with no region code attached.

import os
import json
import urllib.request
import pandas as pd
import geopandas as gpd
import mapclassify
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import folium
from folium import plugins
import branca.colormap as cm
from matplotlib_scalebar.scalebar import ScaleBar
NATIONAL_RATE = 21.3
BREAKS = [10.6, 15, 20, 21.3, 25, 32]
CLASS_COLORS = ["darkgreen", "yellowgreen", "gold", "indianred", "firebrick"]

os.makedirs("data", exist_ok=True)
boundary_path = "data/mar_adm1.geojson"
api_url = "https://www.geoboundaries.org/api/current/gbOpen/MAR/ADM1/"
if not os.path.exists(boundary_path):
    try:
        with urllib.request.urlopen(api_url, timeout=30) as response:
            record = json.load(response)
        print("geoBoundaries API record:", record["boundaryName"], record["boundaryYearRepresented"], record["boundaryLicense"])
        urllib.request.urlretrieve(record["simplifiedGeometryGeoJSON"], boundary_path)
    except Exception as error:
        print("API not reachable, using the GitHub release instead:", error)
        github_url = "https://media.githubusercontent.com/media/wmgeolab/geoBoundaries/main/releaseData/gbOpen/MAR/ADM1/geoBoundaries-MAR-ADM1_simplified.geojson"
        urllib.request.urlretrieve(github_url, boundary_path)
regions = gpd.read_file(boundary_path)
regions = regions[["shapeISO", "shapeName", "geometry"]]
regions = regions.rename(columns={"shapeISO": "iso_code", "shapeName": "boundary_name"})
print("Polygons:", len(regions))
print("CRS:", regions.crs)

Polygons: 12
CRS: EPSG:4326

census = pd.read_csv("data/regions_rgph2024.csv")
capitals = pd.read_csv("data/regional_capitals.csv")
print("Census rows:", len(census))
print(f"Sum of the 12 regional populations: {census['population_2024'].sum():,}")

Census rows: 12
Sum of the 12 regional populations: 36,828,330

The 12 regional populations add up to 36,828,330, the national legal population HCP published, so no region is missing from the census table.

Coordinate reference systems

Two coordinate reference systems are used, each for one job.

• EPSG:4326 (WGS 84 latitude and longitude) is the system the boundary file and the capitals file arrive in. Both joins happen in it, and Folium needs it for the web map.
• ESRI:102022 (Africa Albers Equal Area Conic) preserves area. Every area calculation and the static map use it. Areas computed in EPSG:4326 would be in square degrees, and a square degree near Dakhla (23 degrees north) covers more ground than one near Tangier (36 degrees north).

capital_points = gpd.GeoDataFrame(capitals, geometry=gpd.points_from_xy(capitals["lon"], capitals["lat"]), crs="EPSG:4326")
equal_area = "ESRI:102022"
print("Regions CRS:", regions.crs.to_string())
print("Capitals CRS:", capital_points.crs.to_string())
print("Same CRS before joining:", regions.crs == capital_points.crs)
print("Equal area CRS:", gpd.GeoSeries([], crs=equal_area).crs.name)

Regions CRS: EPSG:4326
Capitals CRS: EPSG:4326
Same CRS before joining: True
Equal area CRS: Africa_Albers_Equal_Area_Conic

Attribute join

The census table and the polygons share the ISO code (MA-01 to MA-12). The join is a left join from the polygons, so every polygon is kept, and validate="one_to_one" stops the notebook if a code appears twice on either side. The two checks after it list any code that exists on only one side.

gdf = regions.merge(census, on="iso_code", how="left", validate="one_to_one")
only_in_boundaries = set(regions["iso_code"]) - set(census["iso_code"])
only_in_census = set(census["iso_code"]) - set(regions["iso_code"])
print("Rows after join:", len(gdf))
print("Codes in boundaries but not in census:", only_in_boundaries or "none")
print("Codes in census but not in boundaries:", only_in_census or "none")
assert len(gdf) == 12
assert gdf["population_2024"].notna().all()
gdf[["iso_code", "boundary_name", "region"]].head(4)

Rows after join: 12
Codes in boundaries but not in census: none
Codes in census but not in boundaries: none

	iso_code	boundary_name	region
0	MA-11	Laâyoune-Sakia El Hamra	Laayoune-Sakia El Hamra
1	MA-01	Tangier-Tetouan-Al Hoceima	Tanger-Tetouan-Al Hoceima
2	MA-02	Oriental	Oriental
3	MA-08	Drâa-Tafilalet	Draa-Tafilalet

The join matched all 12 polygons. The boundary file spells names with accents (Drâa-Tafilalet) and the census file without them (Draa-Tafilalet), which is why the join runs on the ISO code and never on the name.

Normalization

Raw population is a poor variable to map. Large desert regions such as Dakhla-Oued Ed-Dahab would look empty because they are large, and Casablanca-Settat would swamp the scale because it is crowded. Two normalized variables are used instead.

• Unemployment rate is already normalized by HCP: it is the share of the active population aged 15 and over that is out of work and looking. It can be mapped directly.
• Population density divides population by land area in square kilometres, computed in the equal area CRS.

gdf["area_km2"] = gdf.to_crs(equal_area).geometry.area / 1000000
gdf["density_per_km2"] = gdf["population_2024"] / gdf["area_km2"]
gdf["pop_share_pct"] = 100 * gdf["population_2024"] / gdf["population_2024"].sum()
gdf["gap_vs_national"] = gdf["unemployment_rate_2024"] - NATIONAL_RATE
gdf[["region", "population_2024", "area_km2", "density_per_km2"]].round(1).sort_values("density_per_km2", ascending=False)

	region	population_2024	area_km2	density_per_km2
6	Casablanca-Settat	7688967	20297.3	378.8
9	Rabat-Sale-Kenitra	5132639	17630.9	291.1
1	Tanger-Tetouan-Al Hoceima	4030222	16167.6	249.3
7	Marrakech-Safi	4892393	38995.8	125.5
10	Fes-Meknes	4467911	38862.1	115.0
11	Beni Mellal-Khenifra	2525801	27656.0	91.3
4	Souss-Massa	3020431	53638.3	56.3
2	Oriental	2294665	66099.5	34.7
3	Draa-Tafilalet	1655623	86109.2	19.2
5	Guelmim-Oued Noun	448685	44655.7	10.0
0	Laayoune-Sakia El Hamra	451028	144161.0	3.1
8	Dakhla-Oued Ed-Dahab	219965	128508.8	1.7

Classification

With only 12 regions, the choice of class breaks changes the map. The cell below compares two standard schemes with the fixed breaks used here.

rates = gdf["unemployment_rate_2024"]
print("Quantiles (5 classes):", mapclassify.Quantiles(rates, k=5).bins.round(1))
print("Natural breaks (5 classes):", mapclassify.NaturalBreaks(rates, k=5).bins.round(1))
fixed = mapclassify.UserDefined(rates, bins=BREAKS[1:])
print("Fixed breaks used:", BREAKS[1:])
print("Regions per fixed class:", [int(c) for c in fixed.counts])

Quantiles (5 classes): [19.6 20.2 22.9 26.8 31.5]
Natural breaks (5 classes): [10.6 20.8 23.3 26.8 31.5]
Fixed breaks used: [15, 20, 21.3, 25, 32]
Regions per fixed class: [1, 4, 1, 2, 4]

Quantiles force about two regions into each class, so a 0.5 point difference can look as large as a 10 point one. Natural breaks depend on this particular set of 12 values and would move with the next census. The fixed breaks are centred on the national rate of 21.3 percent: the two green classes sit below it, the gold class is the band from 20.0 to 21.3 just under it, and the two red classes sit above it. A reader can tell from the colour alone whether a region is above or below the national figure, which is the question the map asks.

Regional figures

ranked_table = gdf.sort_values("unemployment_rate_2024", ascending=False)
cells = []
cell_colors = []
for index, row in ranked_table.iterrows():
    cells.append([row["region"], f"{row['population_2024']:,.0f}", f"{row['pop_share_pct']:.1f}", f"{row['density_per_km2']:.1f}", f"{row['unemployment_rate_2024']:.1f}", f"{row['gap_vs_national']:+.1f}"])
    shade = "mistyrose" if row["gap_vs_national"] > 0 else "honeydew"
    cell_colors.append([shade] * 6)
headers = ["Region", "Population", "Share of\npopulation (%)", "Inhabitants\nper km²", "Unemployment\n(%)", "Gap vs national\n(points)"]
fig, ax = plt.subplots(figsize=(9, 5))
ax.set_axis_off()
tab = ax.table(cellText=cells, colLabels=headers, cellColours=cell_colors, colColours=["firebrick"] * 6, colWidths=[0.27, 0.13, 0.15, 0.14, 0.15, 0.16], loc="center", cellLoc="center")
tab.auto_set_font_size(False)
tab.set_fontsize(9)
tab.scale(1, 1.6)
for (r, c), cell in tab.get_celld().items():
    cell.set_edgecolor("goldenrod")
    if r == 0:
        cell.set_text_props(color="white", fontweight="bold", fontsize=8)
        cell.set_height(cell.get_height() * 1.6)
    if c == 0 and r > 0:
        cell.set_text_props(ha="left")
plt.show()

Figure 1: Regions ranked by census unemployment rate, RGPH 2024. Red rows are above the national rate of 21.3 percent, green rows below. Density uses area computed in Africa Albers Equal Area. Source: HCP (2024); geoBoundaries (Runfola et al., 2020).

Spatial join

The capitals file has coordinates only. A point in polygon join assigns each capital to the region that contains it. Both layers are in EPSG:4326, checked above, so the join compares like with like. The checks confirm that every capital found a region and that no region received two capitals.

points = gpd.sjoin(capital_points, gdf[["iso_code", "region", "population_2024", "unemployment_rate_2024", "geometry"]], how="left", predicate="within")
print("Capitals matched:", points["region"].notna().sum(), "of", len(points))
print("Distinct regions matched:", points["region"].nunique())
assert points["region"].notna().all()
assert points["region"].is_unique
points[["city", "region", "unemployment_rate_2024"]]

Capitals matched: 12 of 12
Distinct regions matched: 12

	city	region	unemployment_rate_2024
0	Tangier	Tanger-Tetouan-Al Hoceima	19.6
1	Oujda	Oriental	30.4
2	Fes	Fes-Meknes	23.3
3	Rabat	Rabat-Sale-Kenitra	19.8
4	Beni Mellal	Beni Mellal-Khenifra	26.8
5	Casablanca	Casablanca-Settat	18.8
6	Marrakech	Marrakech-Safi	20.8
7	Errachidia	Draa-Tafilalet	22.2
8	Agadir	Souss-Massa	19.7
9	Guelmim	Guelmim-Oued Noun	31.5
10	Laayoune	Laayoune-Sakia El Hamra	26.6
11	Dakhla	Dakhla-Oued Ed-Dahab	10.6

Static map

gdf_eq = gdf.to_crs(equal_area)
points_eq = points.to_crs(equal_area)
legend_labels = ["10.6 to 15.0 (well below)", "15.0 to 20.0 (below)", "20.0 to 21.3 (just below)", "21.3 to 25.0 (above)", "25.0 to 31.5 (well above)"]
fig, ax = plt.subplots(figsize=(8, 9))
gdf_eq.plot(column="unemployment_rate_2024", cmap=ListedColormap(CLASS_COLORS), scheme="UserDefined", classification_kwds={"bins": BREAKS[1:]}, edgecolor="white", linewidth=0.8, legend=True, legend_kwds={"title": "Unemployment rate (%)", "loc": "lower right", "labels": legend_labels, "frameon": False}, ax=ax)
gdf_eq.dissolve().boundary.plot(ax=ax, color="black", linewidth=0.8)
points_eq.plot(ax=ax, color="gold", edgecolor="black", markersize=30, zorder=3)
label_offsets = {"Casablanca": (-50, -6)}
for index, row in points_eq.iterrows():
    offset = label_offsets.get(row["city"], (5, 4))
    ax.annotate(row["city"], xy=(row.geometry.x, row.geometry.y), xytext=offset, textcoords="offset points", fontsize=8, bbox={"boxstyle": "round,pad=0.15", "facecolor": "white", "edgecolor": "none", "alpha": 0.8})
ax.add_artist(ScaleBar(1, units="m", location="lower left", box_alpha=0))
ax.annotate("N", xy=(0.06, 0.95), xytext=(0.06, 0.87), xycoords="axes fraction", ha="center", fontsize=12, fontweight="bold", arrowprops={"arrowstyle": "-|>", "color": "black"})
ax.set_title("Unemployment rate by region, Morocco, 2024 census", fontsize=14, fontweight="bold", loc="left")
ax.set_axis_off()
plt.show()

Figure 2: Census unemployment rate by region, RGPH 2024 (percent of active population aged 15 and over). Fixed breaks are centred on the national rate of 21.3 percent: green below, gold just under, red above. Gold dots are regional capitals. Projection: Africa Albers Equal Area (ESRI:102022). Source: HCP (2024); boundaries from geoBoundaries (Runfola et al., 2020).

The Atlantic core from Tangier to Agadir, where the ports, the automotive plants at Tangier Med and Kenitra and most foreign direct investment sit, stays below the national rate. The interior and the eastern border do worse: the Oriental (30.4 percent) borders a closed Algerian frontier, and Beni Mellal-Khenifra (26.8 percent) depends on agriculture and on remittances. Dakhla-Oued Ed-Dahab (10.6 percent) is the low outlier; it has the smallest population of any region (219,965) and an economy built on fisheries and public investment.

ranked = gdf.sort_values("gap_vs_national")
colors = ["firebrick" if g > 0 else "darkgreen" for g in ranked["gap_vs_national"]]
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.hlines(y=ranked["region"], xmin=0, xmax=ranked["gap_vs_national"], color=colors, linewidth=2)
ax.scatter(ranked["gap_vs_national"], ranked["region"], color=colors, s=50, zorder=3)
ax.axvline(0, color="goldenrod", linewidth=1.5)
ax.set_xlabel("Percentage points above or below the national rate")
ax.spines[["top", "right", "left"]].set_visible(False)
ax.tick_params(axis="y", length=0)
plt.show()

Figure 3: Gap between each region’s census unemployment rate and the national rate of 21.3 percent, in percentage points. Source: HCP (2024).

Six regions sit above the national rate and six below it. The six below include the three most populous, Casablanca-Settat, Rabat-Sale-Kenitra and Marrakech-Safi, so the national average leans toward the coast.

Interactive map

The interactive map has three layers, listed in the layer control at the top right.

• Unemployment rate is on by default and uses the same five classes and colours as the static map, so the two maps read the same way.
• Population density is off by default; turning it on swaps the choropleth. Its scale is in steps because Casablanca-Settat (about 380 per km²) is more than 200 times denser than Dakhla-Oued Ed-Dahab (under 2 per km²).
• Regional capitals are gold markers. Hovering shows the capital and its region’s rate; clicking opens a popup with population and density.

Hovering over any region shows its name, population, rate, gap to the national rate and density. The zoom and fullscreen buttons are at the bottom right, so the two legends at the top left do not cover them, and the minimap at the bottom left keeps Morocco in view when zoomed in. The basemap is Esri’s light grey canvas, which loads without an API key; the CARTO and OpenStreetMap tiles Folium uses by default now refuse requests from a page opened as a local file.

gdf_web = gdf.to_crs("EPSG:4326")
gdf_web["unemployment_label"] = gdf_web["unemployment_rate_2024"].map("{:.1f} %".format)
gdf_web["gap_label"] = gdf_web["gap_vs_national"].map("{:+.1f} points".format)
gdf_web["density_label"] = gdf_web["density_per_km2"].map("{:,.1f} per km²".format)
gdf_web["population_label"] = gdf_web["population_2024"].map("{:,.0f}".format)
gdf_web = gdf_web[["region", "unemployment_rate_2024", "density_per_km2", "unemployment_label", "gap_label", "density_label", "population_label", "geometry"]]
unemp_scale = cm.StepColormap(CLASS_COLORS, index=BREAKS, vmin=BREAKS[0], vmax=BREAKS[-1], caption="Unemployment rate, RGPH 2024 (%)")
density_scale = cm.StepColormap(["honeydew", "palegreen", "mediumseagreen", "seagreen", "darkgreen"], index=[0, 10, 50, 100, 250, 400], vmin=0, vmax=400, caption="Inhabitants per km², 2024")
fields = ["region", "population_label", "unemployment_label", "gap_label", "density_label"]
aliases = ["Region", "Population", "Unemployment", "Gap vs national", "Density"]
esri_url = "https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}"
esri_attr = "Tiles: Esri, HERE, Garmin, OpenStreetMap contributors"
m = folium.Map(tiles=None, zoom_control="bottomright")
folium.TileLayer(tiles=esri_url, attr=esri_attr, name="Light grey basemap", control=False).add_to(m)
folium.GeoJson(gdf_web, name="Unemployment rate", style_function=lambda f: {"fillColor": unemp_scale(f["properties"]["unemployment_rate_2024"]), "color": "white", "weight": 1, "fillOpacity": 0.8}, highlight_function=lambda f: {"weight": 3, "color": "black"}, tooltip=folium.GeoJsonTooltip(fields=fields, aliases=aliases, sticky=True)).add_to(m)
folium.GeoJson(gdf_web, name="Population density", show=False, style_function=lambda f: {"fillColor": density_scale(min(f["properties"]["density_per_km2"], 399)), "color": "white", "weight": 1, "fillOpacity": 0.8}, highlight_function=lambda f: {"weight": 3, "color": "black"}, tooltip=folium.GeoJsonTooltip(fields=fields, aliases=aliases, sticky=True)).add_to(m)
capital_layer = folium.FeatureGroup(name="Regional capitals")
for index, row in points.iterrows():
    popup_text = f"<b>{row['city']}</b><br>Capital of {row['region']}<br>Regional population: {row['population_2024']:,.0f}<br>Unemployment: {row['unemployment_rate_2024']:.1f} %"
    folium.CircleMarker(location=[row["lat"], row["lon"]], radius=6, color="black", weight=1.5, fill=True, fill_color="gold", fill_opacity=1, tooltip=f"{row['city']}: {row['unemployment_rate_2024']:.1f} %", popup=folium.Popup(popup_text, max_width=250)).add_to(capital_layer)
capital_layer.add_to(m)
unemp_scale.add_to(m)
density_scale.add_to(m)
plugins.Fullscreen(position="bottomright").add_to(m)
plugins.MiniMap(tile_layer=folium.TileLayer(tiles=esri_url, attr=esri_attr), toggle_display=True, position="bottomleft").add_to(m)
folium.LayerControl(collapsed=False).add_to(m)
bounds = gdf_web.total_bounds
m.fit_bounds([[bounds[1], bounds[0]], [bounds[3], bounds[2]]])
m

Make this Notebook Trusted to load map: File -> Trust Notebook

Figure 4: Interactive map of census unemployment and population density by region, RGPH 2024, with regional capitals. Source: HCP (2024); geoBoundaries (Runfola et al., 2020).

Switching to the density layer shows that unemployment does not track crowding. Casablanca-Settat is the densest region and has the second lowest rate; the Oriental and Guelmim-Oued Noun are thinly settled and have the two highest.

Reflection and limitations

The maps show a split between the coast and the interior, but the design limits what that split can prove. The first limit is scale. Twelve regions average away most of the variation inside them: Fes-Meknes at 23.3 percent combines the city of Fes with the rural Middle Atlas, and nothing on the map says where inside the region the unemployed live. This is the modifiable areal unit problem. The same census counts aggregated to Morocco’s 75 provinces and prefectures could draw a different pattern, and the regional map would give no warning of it.

Scale also limits what the map says about people. A high regional rate describes a place. Reading from it why any given resident is out of work is the ecological fallacy, and the paragraphs above that link the Oriental’s rate to the closed Algerian border make exactly that kind of claim about a whole region. They are hypotheses the map suggests. Nothing in this notebook tests them, and the policy section that follows should be read the same way.

The indicator itself carries a second set of limits. Census unemployment (21.3 percent nationally) runs well above the labour force survey rate (13.3 percent), and the distance between the two is not the same in every region, so a reader who checks these maps against HCP’s quarterly releases will find different numbers for the same places. Two of the twelve rates are also weaker than the rest. The figures for Marrakech-Safi (20.8) and Draa-Tafilalet (22.2) were read from an HCP chart by rank order and have no second source here, and both sit in classes next to the national rate, so a misreading would change their colour on the map.

Normalization solves one problem and opens another. Dividing population by area removes the size effect, but in the desert regions the result describes almost no one. Laayoune-Sakia El Hamra and Dakhla-Oued Ed-Dahab are mostly empty land with their people concentrated in one or two cities, so a density of 2 or 3 per km² matches nobody’s actual surroundings. The areas themselves come from simplified polygons and do not match HCP’s official surface figures; the Oriental polygon gives about 66,100 km², well under the official surface.

The classification is a choice of its own. Centring the breaks on 21.3 percent answers whether a region is above or below the national rate, and it answers less well how far apart two red regions are. The gap chart carries that second question, which is why the notebook keeps both. The last limit is political. The geoBoundaries file follows Morocco’s administrative division, which includes Laayoune-Sakia El Hamra, Dakhla-Oued Ed-Dahab and part of Guelmim-Oued Noun, while the United Nations lists Western Sahara as a non self governing territory. The map uses these units because they are the units HCP reports the census in.

What the map argues for

• Public job programmes such as Awrach, launched in 2022, allocate places nationally. A regional quota weighted by the census gap (Guelmim-Oued Noun at 10.2 points above the national rate, the Oriental at 9.1) would send money where the shortfall is largest.
• The Oriental’s rate is 11.6 points above Casablanca-Settat’s, and the land border with Algeria has been closed since 1994. The region is the clearest case for testing whether cross border trade, or its absence, explains a regional labour market.
• The next step is the provincial level. HCP publishes census figures for each province, and the same pipeline would run on them with a change of boundary file and the API request for ADM2.

References

geoBoundaries (2026) geoBoundaries API: Morocco ADM1, gbOpen. Available at: https://www.geoboundaries.org/api/current/gbOpen/MAR/ADM1/.

Haut-Commissariat au Plan (2024) Recensement General de la Population et de l’Habitat 2024: principaux resultats. Rabat: HCP, 17 December.

Haut-Commissariat au Plan (2025) Activite, emploi et chomage: resultats annuels 2024. Rabat: HCP, February.

Le Matin (2024) ‘RGPH 2024: le taux de chomage monte a 21,3% en 10 ans’, Le Matin, December. Available at: https://lematin.ma/economie/rgph-2024-le-taux-de-chomage-monte-a-213-en-10-ans/256076.

Runfola, D. et al. (2020) ‘geoBoundaries: a global database of political administrative boundaries’, PLoS ONE, 15(4), e0231866.