This guide has 2 routes. Pick one, or read both. Each route is complete on its own.
Almost every project needs a boundary: the district to study, the state to compare against, the villages to summarise by. The boundary also decides the result. The same district has a different area in different layers, and its neighbours may be named differently. This guide downloads the boundaries of India from three free sources and compares them.
The three sources
| Source | Levels | Strengths | Things to know |
|---|---|---|---|
| geoBoundaries | Country, state, district and more, for every country | An API that returns the source, the year and the licence. A simplified version for maps and a full version for analysis | The detail varies with the country's source data |
| GADM | Country, state, district, sub-district and more | One file per country and level, well known, long history | Free for academic and other non-commercial use, and the terms say you must not redistribute it. Check them if your use is commercial |
| Natural Earth | Country and state, at three scales | Public domain, small files, good for world and country maps | Too coarse for district work |
Which one to pick depends on the job. For district analysis in India, geoBoundaries is a good first choice because the API names its source and year. GADM is the one most often cited in papers. Natural Earth is for maps of a whole country or a continent.
Route A: geoBoundaries
geoBoundaries has an API. You ask for a country and a level, and it answers with information about the layer and the addresses of the files.
You need Python with requests and geopandas.
import requests
import geopandas as gpd
# Step 1: ask the API about the district level (ADM2) of India. ADM0 is the country, ADM1 the state, ADM2 the district.
api = "https://www.geoboundaries.org/api/current/gbOpen/IND/ADM2/"
meta = requests.get(api, timeout=60).json()
if isinstance(meta, list): # some answers come as a list of one item
meta = meta[0]
# Step 2: what is this layer? Read the source, the year and the licence before you use it.
for key in ["boundaryName", "boundaryType", "boundaryYearRepresented", "boundarySource", "boundaryLicense", "admUnitCount"]:
print(f"{key:24s} {meta[key]}")
boundaryName India
boundaryType ADM2
boundaryYearRepresented 2021
boundarySource Pathways Data Pvt. Ltd., lgdirectory.gov.in
boundaryLicense Open Data Commons Open Database License 1.0
admUnitCount 736
The answer is the label that you should put in your report: who made it, for which year, under which licence, and how many units it has. A number of units that does not match what you expect tells you that the data is from a different year.
The API gives two kinds of file. gjDownloadURL is the full detail, about 48 MB for India. simplifiedGeometryGeoJSON has fewer vertices and is about 8 MB. Use the simplified file for maps and for any use where a boundary error of a few tens of metres does not matter. Use the full file for measurements.
# Step 1: read the simplified district file directly from its address. GeoPandas downloads it.
districts = gpd.read_file(meta["simplifiedGeometryGeoJSON"])
print(districts.shape[0], "districts | CRS: EPSG", districts.crs.to_epsg())
print(districts.columns.tolist())
# Step 2: choose one district by name
ludhiana = districts[districts["shapeName"] == "Ludhiana"]
print("found:", len(ludhiana), "| area:", round(ludhiana.to_crs(32643).area.iloc[0] / 1e6, 1), "km2")
735 districts | CRS: EPSG 4326
['shapeName', 'shapeISO', 'shapeID', 'shapeGroup', 'shapeType', 'geometry']
found: 1 | area: 3696.5 km2
The API says 736 units and the simplified file has 735 rows, so even the label and the file do not match exactly. Count what you really got. The names are in the column shapeName. The file has one row for each district and is in latitude and longitude, EPSG:4326. The area was calculated in the UTM zone of the place, because an area in degrees means nothing.
Find the neighbours
A boundary layer also tells you what lies next to a place.
# Step 1: the districts whose outline touches Ludhiana
neighbours = districts[districts.touches(ludhiana.geometry.iloc[0])]
print(len(neighbours), "neighbours:", sorted(neighbours["shapeName"]))
7 neighbours: ['Barnala', 'Fatehgarh Sahib', 'Jalandhar', 'Moga', 'Rupnagar', 'Sangrur', 'Shahid Bhagat Singh Nagar']
States
For a state, ask the API for ADM1 in the same way.
# Step 1: the state level has its own address
meta1 = requests.get("https://www.geoboundaries.org/api/current/gbOpen/IND/ADM1/", timeout=60).json()
meta1 = meta1[0] if isinstance(meta1, list) else meta1
states = gpd.read_file(meta1["simplifiedGeometryGeoJSON"])
punjab = states[states["shapeName"] == "Punjab"]
print(len(states), "states and union territories | Punjab area:", round(punjab.to_crs(32643).area.sum() / 1e6, 1), "km2")
36 states and union territories | Punjab area: 50321.2 km2
Sub-districts
India has a level below the district. geoBoundaries has it as ADM3, and for India it is the sub-district, which is called a tehsil or a taluka in different states. The layer comes from the same source, the Local Government Directory of India, and has about 6,800 units.
# Step 1: ask the API for the sub-district level and read its label first
meta3 = requests.get("https://www.geoboundaries.org/api/current/gbOpen/IND/ADM3/", timeout=60).json()
if isinstance(meta3, list):
meta3 = meta3[0]
for key in ["boundaryType", "boundaryCanonical", "boundaryYearRepresented", "boundarySource", "admUnitCount"]:
print(f"{key:24s} {meta3[key]}")
# Step 2: download the file. It holds the whole country, so it is large, and it takes a while.
subdistricts = gpd.read_file(meta3["gjDownloadURL"])
print(len(subdistricts), "sub-districts | columns:", [c for c in subdistricts.columns if c != "geometry"])
boundaryType ADM3
boundaryCanonical Sub-District
boundaryYearRepresented 2018
boundarySource Pathways Data Pvt. Ltd., lgdirectory.gov.in
admUnitCount 6836
6824 sub-districts | columns: ['shapeName', 'shapeISO', 'shapeID', 'shapeGroup', 'shapeType']
# Step 1: keep the sub-districts whose centre lies inside Ludhiana district. The district file is in the shared data of this section.
ludhiana = gpd.read_file("ludhiana_workflow_data/ludhiana_district.geojson")
centres = subdistricts.to_crs(32643).centroid.to_crs(4326)
inside = subdistricts[centres.within(ludhiana.geometry.iloc[0])].copy()
# Step 2: areas in square kilometres, from the UTM zone of the place
inside["area_km2"] = inside.to_crs(32643).area / 1e6
print(len(inside), "sub-districts in Ludhiana district")
print(inside[["shapeName", "area_km2"]].sort_values("area_km2", ascending=False).round(0).to_string(index=False))
7 sub-districts in Ludhiana district
shapeName area_km2
Jagraon 810.0
Ludhiana (East) 713.0
Samrala 534.0
Ludhiana (West) 529.0
Raikot 460.0
Payal 444.0
Khanna 223.0
The file holds 6,824 sub-districts, a few fewer than the 6,836 that the label counts. The sub-district layer is from 2018 and the district layer used earlier is from 2021, so they come from different years, and where a district was split in between, the two layers will not line up. Here the seven sub-districts add up to about 3,713 square kilometres, close to the 3,697 of the district boundary, and the small difference is where the two sets of edges differ.
Blocks and villages
India has two more levels in geoBoundaries below the sub-district. ADM4 is the community development block, 7,152 units from 2019, and ADM5 is the village, 649,771 units from the 2011 Census. Both are under the same open licence as the other levels. The village file is large: the simplified version is 467 MB and the full one is about 1.1 GB, so this guide uses the simplified one. The block file is 41 MB.
# Step 1: the two levels below the sub-district. Read the label of each before you use it.
for level in ["ADM4", "ADM5"]:
meta = requests.get(f"https://www.geoboundaries.org/api/current/gbOpen/IND/{level}/", timeout=60).json()
if isinstance(meta, list):
meta = meta[0]
print(level, "|", meta["boundaryCanonical"], "|", meta["boundaryYearRepresented"], "|", meta["admUnitCount"], "units |", meta["boundarySource"])
if level == "ADM4":
blocks_url = meta["simplifiedGeometryGeoJSON"]
else:
villages_url = meta["simplifiedGeometryGeoJSON"]
ADM4 | CD Block | 2019 | 7152 units | Pathways Data Pvt. Ltd., lgdirectory.gov.in
ADM5 | village | 2011 | 649771 units | Pathways Data Pvt. Ltd., https//censusindia.gov.in
def fetch(url, name):
# Step 1: download in pieces, and skip the download if the file is already there
path = Path(name)
if not path.exists():
with requests.get(url, stream=True, timeout=300) as reply:
reply.raise_for_status()
with open(path, "wb") as out:
for chunk in reply.iter_content(chunk_size=1 << 22):
out.write(chunk)
return path
blocks_file = fetch(blocks_url, "ind_adm4_simplified.geojson")
villages_file = fetch(villages_url, "ind_adm5_simplified.geojson")
print(blocks_file.name, round(blocks_file.stat().st_size / 1e6), "MB |", villages_file.name, round(villages_file.stat().st_size / 1e6), "MB")
ind_adm4_simplified.geojson 41 MB | ind_adm5_simplified.geojson 467 MB
# Step 1: read only the features whose box meets the district box. This saves memory with a file this large.
ludhiana = gpd.read_file("ludhiana_workflow_data/ludhiana_district.geojson")
centre_of = lambda layer: layer.to_crs(32643).centroid.to_crs(4326)
box = tuple(ludhiana.total_bounds)
blocks = gpd.read_file(blocks_file, bbox=box, engine="pyogrio")
villages = gpd.read_file(villages_file, bbox=box, engine="pyogrio")
# Step 2: keep the units whose centre lies inside the district
blocks = blocks[centre_of(blocks).within(ludhiana.geometry.iloc[0])].copy()
villages = villages[centre_of(villages).within(ludhiana.geometry.iloc[0])].copy()
print(len(blocks), "blocks and", len(villages), "villages with their centre in Ludhiana district")
# Step 3: the size of a village, in hectares
villages["area_ha"] = villages.to_crs(32643).area / 1e4
print("village area: median", round(float(villages["area_ha"].median())), "ha | smallest", round(float(villages["area_ha"].min()), 1), "| largest", round(float(villages["area_ha"].max())), "ha")
print("sum of village areas:", round(float(villages["area_ha"].sum() / 100)), "km2 of", round(float(ludhiana.to_crs(32643).area.iloc[0] / 1e6)), "km2")
print(blocks["shapeName"].tolist())
print(villages["shapeName"].head(5).tolist())
13 blocks and 926 villages with their centre in Ludhiana district
village area: median 295 ha | smallest 16.1 | largest 9955 ha
sum of village areas: 3711 km2 of 3697 km2
['Maloud', 'Ludhiana-1', 'Jagraon', 'Dehlon', 'Machhiwara', 'Sidhwan Bet', 'Doraha', 'Khanna', 'Pakhowal', 'Raikot', 'Samrala', 'Sudhar', 'Ludhiana-2']
['Labangarh', 'Chupki', 'Chak Slhandaud', 'Payal(R)', 'Jaspal Kadar']
Ludhiana district has 13 blocks and 926 villages. The village outlines add up to 3,711 square kilometres, close to the 3,697 of the district boundary, and the median village covers 295 hectares. The two files come from different years, 2019 for blocks and 2011 for villages, and the district boundary is from 2021, so edges do not match exactly, and a village that lay on a district edge was counted in the district that holds its centre. The villages are Census 2011 outlines. Since then, villages have been merged, split and renamed, and some of the names keep marks copied from the Census. Use the layer for the village pattern of 2011, and check a few villages that you know on a map before you rely on it.
The Census of India and the Local Government Directory also publish tables of village names and codes, which are not shapes. You can join such a table to this layer by name, and the join needs the care described in How to turn a table into points and join it to polygons, because spellings differ.
Route B: GADM
GADM serves one file for each country and level. The district level for India is the second level, called level 2. The JSON version is a small zip file.
import io
import zipfile
# Step 1: download the zip of the district level for India, and open the one file inside it
reply = requests.get("https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_IND_2.json.zip", timeout=120)
archive = zipfile.ZipFile(io.BytesIO(reply.content))
print("downloaded", round(len(reply.content) / 1e6, 1), "MB |", archive.namelist())
# Step 2: read it. GADM names the state NAME_1 and the district NAME_2.
gadm = gpd.read_file(archive.open(archive.namelist()[0]))
print(gadm.shape[0], "districts |", gadm.columns.tolist()[:8])
# Step 3: the same district
gadm_ludhiana = gadm[(gadm["NAME_1"] == "Punjab") & (gadm["NAME_2"] == "Ludhiana")]
print("GADM area:", round(gadm_ludhiana.to_crs(32643).area.iloc[0] / 1e6, 1), "km2")
downloaded 1.5 MB | ['gadm41_IND_2.json']
676 districts | ['GID_2', 'GID_0', 'COUNTRY', 'GID_1', 'NAME_1', 'NL_NAME_1', 'NAME_2', 'VARNAME_2']
GADM area: 3566.2 km2
GADM names the state in one column and the district in another, so a name like Aurangabad, which exists in two states, needs both to select the right one.
Route C: Natural Earth
Natural Earth has states at three scales. The 1:10 million version is the most detailed.
# Step 1: read the 10 m scale states and provinces of the world straight from the zip file
world_states = gpd.read_file("zip+https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip")
print(world_states.shape[0], "states and provinces in the world")
# Step 2: pick Punjab, India. The country is in the column admin and the state in name.
ne_punjab = world_states[(world_states["admin"] == "India") & (world_states["name"] == "Punjab")]
print("Natural Earth area of Punjab:", round(ne_punjab.to_crs(32643).area.sum() / 1e6, 1), "km2")
4596 states and provinces in the world
Natural Earth area of Punjab: 50344.7 km2
Natural Earth does not have districts. It is a map for looking at, and its borders are drawn at a scale of 1:10 million, so use it for context and for national and state maps and not for measurements.
How much do the sources disagree?
import pandas as pd
# Step 1: the area of Ludhiana district in the two district sources, and of Punjab in the two that have states
table = pd.DataFrame({
"layer": ["Ludhiana district", "Ludhiana district", "Punjab state", "Punjab state"],
"source": ["geoBoundaries", "GADM", "geoBoundaries", "Natural Earth"],
"area_km2": [
ludhiana.to_crs(32643).area.iloc[0] / 1e6, gadm_ludhiana.to_crs(32643).area.iloc[0] / 1e6,
punjab.to_crs(32643).area.sum() / 1e6, ne_punjab.to_crs(32643).area.sum() / 1e6,
],
}).round(1)
print(table.to_string(index=False))
layer source area_km2
Ludhiana district geoBoundaries 3696.5
Ludhiana district GADM 3566.2
Punjab state geoBoundaries 50321.2
Punjab state Natural Earth 50344.7
The two district layers differ by about 130 km², roughly 3.5 percent, for the same district. The two state layers differ by less than 0.1 percent. The district difference is large enough to matter for an area figure or a density. It comes from different source data and different years. geoBoundaries describes its India district layer as 736 units from 2021, and the GADM file has 676 districts, which suggests that districts were created or reorganised between the two. A boundary layer is a snapshot of a map that keeps changing.
The practical rules follow from this:
- Choose one source for a project, and use it for every layer. Do not mix a GADM boundary with a geoBoundaries one.
- Write the source, the year and the licence in your report.
- Quote an area from your layer as the area of that layer, and not as the official area of the district.
- When the names have to match another table, such as census data, join on a code and not on the name. Spellings differ.
Save a boundary
# Step 1: write the district to a GeoPackage in its own CRS, and the same in a UTM zone for measurements
ludhiana.to_file("ludhiana_district.gpkg", layer="district_4326", driver="GPKG")
ludhiana.to_crs(32643).to_file("ludhiana_district.gpkg", layer="district_utm43n", driver="GPKG")
# Step 2: read the layers back to check
print(gpd.list_layers("ludhiana_district.gpkg"))
name geometry_type
0 district_4326 Polygon
1 district_utm43n Polygon
Before you publish
- Licence. geoBoundaries gbOpen is under the Open Data Commons Open Database Licence. Credit the project and its source. GADM has its own terms of use and does not permit redistributing the data. Natural Earth is public domain.
- Borders. The international and disputed borders of a country are drawn differently by different sources, and India has rules about how its borders may appear on published maps. For a map that you publish, check the current rules, and use a boundary that follows them. The free layers above are for analysis and for study use.
- Detail. A simplified boundary is fine for a map at the scale of a district. It is not for a parcel or a village.
Common mistakes
- Measuring an area in degrees. Reproject to a UTM zone first.
- Joining a table to the boundary by district name and losing the districts that are spelled differently.
- Mixing boundaries from two sources in one map, and getting slivers and gaps where they disagree.
- Using the boundary of a year that does not match your data, in a country where districts change.
- Treating an area from a boundary layer as the official figure.
- Forgetting the licence and the credit line.