[Python] Why does my shapefile color a different polygon than intended?

I have risk data which I would like to color on a map according to the risk level. I read a shapefile and a csv data file which I then merge. This works very well working with adm1 shapefile.

When I run the same script with adm2 shapefile, the results are completely weird: The polygons colored are way far from the polygon with data. I have the attached map 3 examples of the colored polygons and in black dot the location of the data.

I will appreciate if anyone can give a clue of what is going on, and if possible how to resolve this.

Python script used to produce the plot:

#!/home/zmumba/anaconda3/bin/python

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from mpl_toolkits.basemap import Basemap # optional ?

map_df = gpd.read_file("/home/zmumba/DA/Dekad_Progs/Shapefiles/Lesotho/geoBoundaries-LSO-ADM2-all/geoBoundaries-LSO-ADM2.shx")
risks_df=pd.read_csv("/home/zmumba/DA/Dekad_Progs/Output/H1dRrisks.csv")

map_df["risk"] = map_df.merge(risks_df, left_on="shapeName", right_on="District")["risk"]

colors = {1: "green", 2: "yellow", 3: "orange", 4: "red"} # or a list
labels = {1: "no risk", 2: "low risk", 3: "medium risk", 4: "high risk"}
catego = map_df["risk"].astype(str).str.cat(map_df["risk"].map(labels), sep="- ")

fig, ax = plt.subplots(figsize=(5, 5))
plt.title(f'Risk of Heavy 24hr Rain: 20-24Nov', y=1.04)

map_df.plot(
column=catego,
categorical=True,
edgecolor="k",
linewidths=0.8,
alpha=0.7,
cmap=ListedColormap([c for r,c in colors.items() if r in map_df["risk"].unique()]),
legend=True,
legend_kwds={
"title": "Risk Level",
"shadow": True,
"loc": "lower right",
"fontsize": 10,
},
ax=ax,
)
ax.set_axis_off()
plt.savefig('H1dRriskmap.png', dpi=300)

plt.show()


Just to add some clarification to my problem, The attached image shows the polygon numbers from the shapefile. Coloring 42 colors 41, coloring 43 colors 42. The coordinates of 41, 42 and 43 are

-29.3892289999999, 28.3056183629316 -> 41
-29.2870792999999, 29.2715247780569 -> 42
-29.2255776170000, 27.6546474532047 -> 43


These coordinates are the centroids of the polygons taken from the shapefile itself. So it cannot be a problem of coordinates. Could there be something wrong in the python code? The code reads from a file which is in the format:

"District","risk"
"name1",1
"name2",1
...
"name78",1


where 1 can b2 2, 3, or 4 depending on the risk level. I would be glad to try doing the same thing in R, but the problem with R is, it will say "no package sf" and so on.

Continue reading...