dc-l2-cairo-invalid-dedup
Map (reference ← swipe → agent)
0:00
cairo_parcels_legacy.geojson is a parcel snapshot stitched together from three legacy provincial systems before unification, and I need to turn it into the foundation for our new central registry. Can you clean it up so we end up with one canonical record per parcel? When duplicates exist, keep the record with the earliest record_seq, and discard any polygon fragments smaller than 1 m².
Please write the result to parcels_canonical.geoparquet including these columns: parcel_id (the join key), parcel_class, district, and area_m2. Every feature must be a MultiPolygon, and there should be no null or empty geometries, no invalid geometries, no duplicate geometries, and no artifact slivers.
0:00
System info
{
"adapter": "openrouter",
"model": "deepseek/deepseek-v4-flash",
"base_url": "https://openrouter.ai/api/v1"
}
0:07
Let me start by reading the input file to understand its structure and content.
0:07
Read
/work/cairo_parcels_legacy.geojson
1 {
2 "type": "FeatureCollection",
3 "name": "cairo_parcels_legacy",
4 "crs": {
5 "type": "name",
6 "properties": {
7 "name": "urn:ogc:def:crs:EPSG::22992"
8 }
9 },
10 "features": [
11 {
12 "type": "Feature",
13 "geometry": {
14 "type": "Polygon",
15 "coordinates": [
16 [
17 [
18 640000.0,
19 815000.0
20 ],
21 [
22 640030.0,
23 815000.0
24 ],
25 [
26 640030.0,
27 815040.0
28 ],
29 [
30 640000.0,
31 815040.0
32 ],
33 [
34 640000.0,
35 815000.0
36 ]
37 ]
38 ]
39 },
40 "properties": {
41 "parcel_id": 1,
42 "record_seq": 1,
43 "parcel_class": "commercial",
44 "district": "Giza-East",
45 "area_m2": 1200.0
46 }
47 },
48 {
49 "type": "Feature",
50 "geometry": {
51 "type": "Polygon",
52 "coordinates": [
53 [
54 [
55 640050.0,
56 815000.0
57 ],
58 [
59 640080.0,
60 815000.0
<system-reminder>Showed lines 1-60 of 11522. Re-call Read with offset=61 to continue, or raise limit.</system-reminder>
0:13
Let me understand the data: 1. It's a GeoJSON file with EPSG:22992 CRS (Egyptian Red Belt / Modified, metric coordinates) 2. It has features with properties: parcel_id, record_seq, parcel_class, district, area_m2 3. I need to: - Keep one record per parcel_id (the one with earliest record_seq) - Remove polygon fragments < 1 m² - Ensure geometries are MultiPolygon (not Polygon) - No null/empty/invalid geometries - No duplicate geometries - No slivers - Output: parcels_canonical.geoparquet with columns: parcel_id, parcel_class, district, area_m2 Let me first explore the data to understand its structure and issues.
0:13
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import json
with open('cairo_parcels_legacy.geojson') as f:
data = json.load(f)
print('Total features:', len(data['features']))
props = [f['properties'] for f in data['features']]
keys = props[0].keys()
print('Property keys:', list(keys))
# Check geometry types
geom_types = {}
for f in data['features']:
gt = f['geometry']['type'] if f['geometry'] else 'None'
geom_types[gt] = geom_types.get(gt, 0) + 1
print('Geometry types:', geom_types)
# Parcel IDs distribution
from collections import Counter
pids = Counter(p['parcel_id'] for p in props)
print('Unique parcel_ids:', len(pids))
print('Parcel IDs with >1 record:', sum(1 for v in pids.values() if v > 1))
# record_seq range
seqs = [p['record_seq'] for p in props]
print('record_seq range:', min(seqs), '-', max(seqs))
"
Total features: 290
Property keys: ['parcel_id', 'record_seq', 'parcel_class', 'district', 'area_m2']
Geometry types: {'Polygon': 260, 'MultiPolygon': 30}
Unique parcel_ids: 290
Parcel IDs with >1 record: 0
record_seq range: 1 - 290
0:19
Interesting - all 290 parcel_ids are unique, so there are no duplicates by parcel_id. But the user says "when duplicates exist, keep the record with the earliest record_seq". Let me look more carefully - perhaps there are parcels with same geometry but different parcel_id? Or the user's wording is just saying "handle duplicates as a general rule". Let me check for duplicate geometries and look at the data more carefully.
0:19
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import json
from shapely import wkt, from_geojson, to_wkt
with open('cairo_parcels_legacy.geojson') as f:
data = json.load(f)
# Check for duplicate geometries
geoms = []
for feat in data['features']:
g = from_geojson(json.dumps(feat['geometry']))
geoms.append(g)
# Check for duplicates
from collections import Counter
wkt_list = [to_wkt(g) for g in geoms]
dup_count = Counter(wkt_list)
print('Duplicate geometries:', sum(1 for v in dup_count.values() if v > 1))
for w, c in dup_count.items():
if c > 1:
print(f' Count {c}: {w[:100]}...')
# Check null/empty geometries
null_or_empty = sum(1 for g in geoms if g is None or g.is_empty)
print('Null/empty geometries:', null_or_empty)
# Check invalid geometries
invalid = sum(1 for g in geoms if g is not None and not g.is_valid)
print('Invalid geometries:', invalid)
# Check area distribution
areas = [g.area for g in geoms if g is not None]
print('Min area:', min(areas))
print('Max area:', max(areas))
small = sum(1 for a in areas if a < 1.0)
print('Areas < 1 m²:', small)
# Check CRS
print('CRS:', data.get('crs'))
" 2>&1
Duplicate geometries: 50
Count 2: POLYGON ((640000 815000, 640030 815000, 640030 815040, 640000 815040, 640000 815000))...
Count 2: POLYGON ((640050 815000, 640080 815000, 640080 815040, 640050 815040, 640050 815000))...
Count 2: POLYGON ((640100 815000, 640130 815000, 640130 815040, 640100 815040, 640100 815000))...
Count 2: POLYGON ((640150 815000, 640180 815000, 640180 815040, 640150 815040, 640150 815000))...
Count 2: POLYGON ((640200 815000, 640230 815000, 640230 815040, 640200 815040, 640200 815000))...
Count 2: POLYGON ((640250 815000, 640280 815000, 640280 815040, 640250 815040, 640250 815000))...
Count 2: POLYGON ((640350 815000, 640380 815000, 640380 815040, 640350 815040, 640350 815000))...
Count 2: POLYGON ((640400 815000, 640430 815000, 640430 815040, 640400 815040, 640400 815000))...
Count 2: POLYGON ((640450 815000, 640480 815000, 640480 815040, 640450 815040, 640450 815000))...
Count 2: POLYGON ((640550 815000, 640580 815000, 640580 815040, 640550 815040, 640550 815000))...
Count 2: POLYGON ((640600 815000, 640630 815000, 640630 815040, 640600 815040, 640600 815000))...
Count 2: POLYGON ((640700 815000, 640730 815000, 640730 815040, 640700 815040, 640700 815000))...
Count 2: POLYGON ((640000 815060, 640030 815060, 640030 815100, 640000 815100, 640000 815060))...
Count 2: POLYGON ((640100 815060, 640130 815060, 640130 815100, 640100 815100, 640100 815060))...
Count 2: POLYGON ((640200 815060, 640230 815060, 640230 815100, 640200 815100, 640200 815060))...
Count 2: POLYGON ((640250 815060, 640280 815060, 640280 815100, 640250 815100, 640250 815060))...
Count 2: POLYGON ((640300 815060, 640330 815060, 640330 815100, 640300 815100, 640300 815060))...
Count 2: POLYGON ((640400 815060, 640430 815060, 640430 815100, 640400 815100, 640400 815060))...
Count 2: POLYGON ((640450 815060, 640480 815060, 640480 815100, 640450 815100, 640450 815060))...
Count 2: POLYGON ((640500 815060, 640530 815060, 640530 815100, 640500 815100, 640500 815060))...
Count 2: POLYGON ((640550 815060, 640580 815060, 640580 815100, 640550 815100, 640550 815060))...
Count 2: POLYGON ((640600 815060, 640630 815060, 640630 815100, 640600 815100, 640600 815060))...
Count 2: POLYGON ((640700 815060, 640730 815060, 640730 815100, 640700 815100, 640700 815060))...
Count 2: POLYGON ((640000 815120, 640030 815120, 640030 815160, 640000 815160, 640000 815120))...
Count 2: POLYGON ((640100 815120, 640130 815120, 640130 815160, 640100 815160, 640100 815120))...
Count 2: POLYGON ((640150 815120, 640180 815120, 640180 815160, 640150 815160, 640150 815120))...
Count 2: POLYGON ((640250 815120, 640280 815120, 640280 815160, 640250 815160, 640250 815120))...
Count 2: POLYGON ((640350 815120, 640380 815120, 640380 815160, 640350 815160, 640350 815120))...
Count 2: POLYGON ((640400 815120, 640430 815120, 640430 815160, 640400 815160, 640400 815120))...
Count 2: POLYGON ((640450 815120, 640480 815120, 640480 815160, 640450 815160, 640450 815120))...
Count 2: POLYGON ((640550 815120, 640580 815120, 640580 815160, 640550 815160, 640550 815120))...
Count 2: POLYGON ((640600 815120, 640630 815120, 640630 815160, 640600 815160, 640600 815120))...
Count 2: POLYGON ((640650 815120, 640680 815120, 640680 815160, 640650 815160, 640650 815120))...
Count 2: POLYGON ((640700 815120, 640730 815120, 640730 815160, 640700 815160, 640700 815120))...
Count 2: POLYGON ((640000 815180, 640030 815180, 640030 815220, 640000 815220, 640000 815180))...
Count 2: POLYGON ((640100 815180, 640130 815180, 640130 815220, 640100 815220, 640100 815180))...
Count 2: POLYGON ((640150 815180, 640180 815180, 640180 815220, 640150 815220, 640150 815180))...
Count 2: POLYGON ((640250 815180, 640280 815180, 640280 815220, 640250 815220, 640250 815180))...
Count 2: POLYGON ((640300 815180, 640330 815180, 640330 815220, 640300 815220, 640300 815180))...
Count 2: POLYGON ((640400 815180, 640430 815180, 640430 815220, 640400 815220, 640400 815180))...
Count 2: POLYGON ((640500 815180, 640530 815180, 640530 815220, 640500 815220, 640500 815180))...
Count 2: POLYGON ((640550 815180, 640580 815180, 640580 815220, 640550 815220, 640550 815180))...
Count 2: POLYGON ((640650 815180, 640680 815180, 640680 815220, 640650 815220, 640650 815180))...
Count 2: POLYGON ((640700 815180, 640730 815180, 640730 815220, 640700 815220, 640700 815180))...
Count 2: POLYGON ((640000 815240, 640030 815240, 640030 815280, 640000 815280, 640000 815240))...
Count 2: POLYGON ((640050 815240, 640080 815240, 640080 815280, 640050 815280, 640050 815240))...
Count 2: POLYGON ((640100 815240, 640130 815240, 640130 815280, 640100 815280, 640100 815240))...
Count 2: POLYGON ((640150 815240, 640180 815240, 640180 815280, 640150 815280, 640150 815240))...
Count 2: POLYGON ((640250 815240, 640280 815240, 640280 815280, 640250 815280, 640250 815240))...
Count 2: POLYGON ((640300 815240, 640330 815240, 640330 815280, 640300 815280, 640300 815240))...
Null/empty geometries: 0
Invalid geometries: 20
Min area: 0.0
Max area: 1280.0
Areas < 1 m²: 50
CRS: {'type': 'name', 'properties': {'name': 'urn:ogc:def:crs:EPSG::22992'}}
0:27
Let me understand the data: 1. 290 features, but parcel_ids are all unique (290 unique IDs). So no duplicate parcel_ids. 2. But there are 50 duplicate geometries (25 pairs, each pair has same geometry but different parcel_id? or same parcel_id?). Wait, "Parcel IDs with >1 record: 0" - so each parcel_id appears exactly once. But there are duplicate geometries - meaning two different parcel_ids can have the same geometry? That seems odd but it's what the data shows. Actually wait - "Parcel IDs with >1 record: 0" means no parcel_id repeats. So each parcel_id is unique. But there are 50 duplicate geometries. This means 50 pairs of features with different parcel_ids but identical geometries. Let me re-examine. The user says "stitched together from three legacy provincial systems before unification" and "When duplicates exist, keep the record with the earliest record_seq." But there are no duplicate parcel_ids. However, there ARE duplicate geometries. So maybe the "duplicates" the user refers to are records that represent the same physical parcel (same geometry) but from different legacy systems, thus having different parcel_ids and record_seq values. So the logic would be: - Group by geometry (since same geometry = same physical parcel) - Within each group, keep the record with the earliest (smallest) record_seq - Discard polygon fragments smaller than 1 m² (area < 1) Also: 20 invalid geometries, 50 geometries with area < 1 m². Let me check which features have area < 1 and which are invalid to understand the overlap. Actually, let me look at the data more carefully. 50 geometries with area 0.0 and 50 duplicate geometries. Let me see if there are small-slop polygons too. Let me write a comprehensive solution.
0:27
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import json
from shapely import from_geojson, to_wkt, is_valid, area
with open('cairo_parcels_legacy.geojson') as f:
data = json.load(f)
# Check which features have 0 area
for i, feat in enumerate(data['features']):
g = from_geojson(json.dumps(feat['geometry']))
a = area(g)
if a == 0.0:
print(f' [{i}] parcel_id={feat[\"properties\"][\"parcel_id\"]}, seq={feat[\"properties\"][\"record_seq\"]}, area={a}, geom={to_wkt(g)[:80]}')
print('---')
# Check invalid
for i, feat in enumerate(data['features']):
g = from_geojson(json.dumps(feat['geometry']))
if not is_valid(g):
print(f' [{i}] parcel_id={feat[\"properties\"][\"parcel_id\"]}, valid={is_valid(g)}, area={area(g)}, geom={to_wkt(g)[:80]}')
" 2>&1
[6] parcel_id=7, seq=7, area=0.0, geom=POLYGON ((640300 815000, 640330 815040, 640300 815040, 640330 815000, 640300 815 [16] parcel_id=17, seq=17, area=0.0, geom=POLYGON ((640050 815060, 640080 815100, 640050 815100, 640080 815060, 640050 815 [22] parcel_id=23, seq=23, area=0.0, geom=POLYGON ((640350 815060, 640380 815100, 640350 815100, 640380 815060, 640350 815 [34] parcel_id=35, seq=35, area=0.0, geom=POLYGON ((640200 815120, 640230 815160, 640200 815160, 640230 815120, 640200 815 [40] parcel_id=41, seq=41, area=0.0, geom=POLYGON ((640500 815120, 640530 815160, 640500 815160, 640530 815120, 640500 815 [52] parcel_id=53, seq=53, area=0.0, geom=POLYGON ((640350 815180, 640380 815220, 640350 815220, 640380 815180, 640350 815 [57] parcel_id=58, seq=58, area=0.0, geom=POLYGON ((640600 815180, 640630 815220, 640600 815220, 640630 815180, 640600 815 [70] parcel_id=71, seq=71, area=0.0, geom=POLYGON ((640500 815240, 640530 815280, 640500 815280, 640530 815240, 640500 815 [75] parcel_id=76, seq=76, area=0.0, geom=POLYGON ((640000 815300, 640030 815340, 640000 815340, 640030 815300, 640000 815 [88] parcel_id=89, seq=89, area=0.0, geom=POLYGON ((640650 815300, 640680 815340, 640650 815340, 640680 815300, 640650 815 [93] parcel_id=94, seq=94, area=0.0, geom=POLYGON ((640150 815360, 640180 815400, 640150 815400, 640180 815360, 640150 815 [106] parcel_id=107, seq=107, area=0.0, geom=POLYGON ((640050 815420, 640080 815460, 640050 815460, 640080 815420, 640050 815 [111] parcel_id=112, seq=112, area=0.0, geom=POLYGON ((640300 815420, 640330 815460, 640300 815460, 640330 815420, 640300 815 [124] parcel_id=125, seq=125, area=0.0, geom=POLYGON ((640200 815480, 640230 815520, 640200 815520, 640230 815480, 640200 815 [129] parcel_id=130, seq=130, area=0.0, geom=POLYGON ((640450 815480, 640480 815520, 640450 815520, 640480 815480, 640450 815 [142] parcel_id=143, seq=143, area=0.0, geom=POLYGON ((640350 815540, 640380 815580, 640350 815580, 640380 815540, 640350 815 [147] parcel_id=148, seq=148, area=0.0, geom=POLYGON ((640600 815540, 640630 815580, 640600 815580, 640630 815540, 640600 815 [160] parcel_id=161, seq=161, area=0.0, geom=POLYGON ((640500 815600, 640530 815640, 640500 815640, 640530 815600, 640500 815 [165] parcel_id=166, seq=166, area=0.0, geom=POLYGON ((640000 815660, 640030 815700, 640000 815700, 640030 815660, 640000 815 [178] parcel_id=179, seq=179, area=0.0, geom=POLYGON ((640650 815660, 640680 815700, 640650 815700, 640680 815660, 640650 815 --- [6] parcel_id=7, valid=False, area=0.0, geom=POLYGON ((640300 815000, 640330 815040, 640300 815040, 640330 815000, 640300 815 [16] parcel_id=17, valid=False, area=0.0, geom=POLYGON ((640050 815060, 640080 815100, 640050 815100, 640080 815060, 640050 815 [22] parcel_id=23, valid=False, area=0.0, geom=POLYGON ((640350 815060, 640380 815100, 640350 815100, 640380 815060, 640350 815 [34] parcel_id=35, valid=False, area=0.0, geom=POLYGON ((640200 815120, 640230 815160, 640200 815160, 640230 815120, 640200 815 [40] parcel_id=41, valid=False, area=0.0, geom=POLYGON ((640500 815120, 640530 815160, 640500 815160, 640530 815120, 640500 815 [52] parcel_id=53, valid=False, area=0.0, geom=POLYGON ((640350 815180, 640380 815220, 640350 815220, 640380 815180, 640350 815 [57] parcel_id=58, valid=False, area=0.0, geom=POLYGON ((640600 815180, 640630 815220, 640600 815220, 640630 815180, 640600 815 [70] parcel_id=71, valid=False, area=0.0, geom=POLYGON ((640500 815240, 640530 815280, 640500 815280, 640530 815240, 640500 815 [75] parcel_id=76, valid=False, area=0.0, geom=POLYGON ((640000 815300, 640030 815340, 640000 815340, 640030 815300, 640000 815 [88] parcel_id=89, valid=False, area=0.0, geom=POLYGON ((640650 815300, 640680 815340, 640650 815340, 640680 815300, 640650 815 [93] parcel_id=94, valid=False, area=0.0, geom=POLYGON ((640150 815360, 640180 815400, 640150 815400, 640180 815360, 640150 815 [106] parcel_id=107, valid=False, area=0.0, geom=POLYGON ((640050 815420, 640080 815460, 640050 815460, 640080 815420, 640050 815 [111] parcel_id=112, valid=False, area=0.0, geom=POLYGON ((640300 815420, 640330 815460, 640300 815460, 640330 815420, 640300 815 [124] parcel_id=125, valid=False, area=0.0, geom=POLYGON ((640200 815480, 640230 815520, 640200 815520, 640230 815480, 640200 815 [129] parcel_id=130, valid=False, area=0.0, geom=POLYGON ((640450 815480, 640480 815520, 640450 815520, 640480 815480, 640450 815 [142] parcel_id=143, valid=False, area=0.0, geom=POLYGON ((640350 815540, 640380 815580, 640350 815580, 640380 815540, 640350 815 [147] parcel_id=148, valid=False, area=0.0, geom=POLYGON ((640600 815540, 640630 815580, 640600 815580, 640630 815540, 640600 815 [160] parcel_id=161, valid=False, area=0.0, geom=POLYGON ((640500 815600, 640530 815640, 640500 815640, 640530 815600, 640500 815 [165] parcel_id=166, valid=False, area=0.0, geom=POLYGON ((640000 815660, 640030 815700, 640000 815700, 640030 815660, 640000 815 [178] parcel_id=179, valid=False, area=0.0, geom=POLYGON ((640650 815660, 640680 815700, 640650 815700, 640680 815660, 640650 815
0:48
Let me understand the data: 1. **290 features** total - 260 Polygons, 30 MultiPolygons 2. **290 unique parcel_ids** - so parcel_id is already unique... but there are duplicate geometries (50 duplicates → 25 pairs) 3. But wait - the user said "When duplicates exist, keep the record with the earliest record_seq" - so duplicates are identified by the same parcel_id? No, parcel_id is already unique (290 unique). Let me check if there are duplicate parcel_ids with different record_seq values... Actually re-reading: "290 unique parcel_ids" means all 290 parcel_ids are distinct, so there are no duplicate parcel_ids. But there are 50 duplicate geometries (25 pairs). So perhaps duplicate means "same geometry" regardless of parcel_id? Let me look more carefully. The duplicate geometries are pairs - same geometry appearing twice with different parcel_id. The user says "When duplicates exist, keep the record with the earliest record_seq". So among records with identical geometry, keep the one with smallest record_seq. Also: - 50 geometries with area < 1 m² (including the 20 invalid ones which are zero-area bowties) - 0 null/empty geometries - All parcel_ids are unique Let me think about the approach: 1. Read the GeoJSON 2. For each feature, compute the geometry's area to filter out < 1 m² 3. For duplicate geometries, keep only the one with smallest record_seq 4. Ensure all geometries are MultiPolygon 5. Ensure all geometries are valid (fix invalid ones, or filter them out since they're zero-area bowties) 6. Write to GeoParquet with columns: parcel_id, parcel_class, district, area_m2 Actually, the zero-area invalid geometries are also < 1 m², so they'd be filtered out by the area filter anyway. Let me write the script carefully. Wait - there are 50 features with area < 1 m²? Let me check more carefully... The output said "Areas < 1 m²: 50" - but the zero-area ones are 20 (all the invalid bowtie ones). Where do the other 30 come from? Actually looking again: "Areas < 1 m²: 50" - so 50 features have area < 1.0. The 20 invalid ones have area 0.0, and there are 30 more small ones. Those 30 are probably the duplicates that are tiny fragments. Wait... the 50 duplicate geometries (25 pairs) - the duplicates have the same geometry. But the area check is against individual geometries. So if there are 25 pairs (50 features total) with duplicate geometries, and their area might be fine... Let me re-examine. The output said "Duplicate geometries: 50" - meaning 50 features have geometries that appear at least twice. Since each appears exactly twice (Count 2), there are 25 pairs = 50 features involved. And "Areas < 1 m²: 50" - 50 features have area less than 1. These could be different features than the duplicates. Let me check by looking at a few features more carefully to understand the structure. Let me just write a comprehensive script that handles everything. Plan: 1. Load the GeoJSON 2. Build a GeoDataFrame 3. Group by geometry (WKT) - for duplicates within same geometry, keep min record_seq 4. Filter out geometries with area < 1 m² 5. Convert all geometries to MultiPolygon 6. Make geometries valid (buffer(0)) 7. Handle any remaining duplicates (after area filter, there might be geometrically identical features) 8. Output to GeoParquet with the required columns Let me also compute area_m2 from the actual geometry rather than using the stored one. Actually wait - the CRS is EPSG:22992 which uses meters, so area from shapely should be in m² directly. Let me verify. EPSG:22992 = Egypt Red Belt / TM. Yes, units are meters. Let me write the script. Actually, I need to check if there are really no duplicate parcel_ids with different record_seq values. The user says "When duplicates exist, keep the record with the earliest record_seq" - suggesting duplicates exist by parcel_id. But from my analysis all parcel_ids are unique... Hmm, maybe "duplicates" refers to duplicate geometries (same shape/shared vertices) rather than duplicate parcel_ids. Or maybe it's both. Let me just handle duplicate geometries as described. Actually, let me re-read: "When duplicates exist, keep the record with the earliest record_seq, and discard any polygon fragments smaller than 1 m²." "Duplicates" could mean features with the same parcel_id (which don't exist), or duplicate geometries. Since all parcel_ids are unique but there are duplicate geometries, I'll interpret "duplicates" as duplicate geometries. Let me write a robust script. Also important: the user says "no duplicate geometries" in the output. So after filtering by area and deduplicating by geometry, we need to make sure no two features have the same geometry. Let me write the script now.
0:48
0:48
Write
/work/solve.py
content · 141 lines · py
import json
import warnings
import geopandas as gpd
import pandas as pd
from shapely import (
from_geojson, to_wkt, area, is_valid, make_valid,
MultiPolygon, Polygon
)
# ---------------------------------------------------------------------------
# 1. Load source
# ---------------------------------------------------------------------------
with open("/work/cairo_parcels_legacy.geojson") as f:
src = json.load(f)
features = src["features"]
print(f"Total input features: {len(features)}")
# ---------------------------------------------------------------------------
# 2. Build a list of dicts for GeoDataFrame construction
# ---------------------------------------------------------------------------
rows = []
for i, feat in enumerate(features):
props = feat["properties"]
geom_json = feat["geometry"]
if geom_json is None:
print(f" SKIP [{i}] null geometry – parcel_id={props['parcel_id']}")
continue
geom = from_geojson(json.dumps(geom_json))
if geom.is_empty:
print(f" SKIP [{i}] empty geometry – parcel_id={props['parcel_id']}")
continue
rows.append({
"parcel_id": props["parcel_id"],
"record_seq": props["record_seq"],
"parcel_class": props["parcel_class"],
"district": props["district"],
"geometry": geom,
})
print(f"After null/empty removal: {len(rows)}")
# ---------------------------------------------------------------------------
# 3. Build GeoDataFrame
# ---------------------------------------------------------------------------
gdf = gpd.GeoDataFrame(rows, geometry="geometry", crs="EPSG:22992")
print(f"CRS: {gdf.crs}")
# ---------------------------------------------------------------------------
# 4. Filter out polygon fragments < 1 m² (before dedup, so the tiny ones
# don't compete with valid duplicates)
# ---------------------------------------------------------------------------
gdf["_area"] = gdf.geometry.area
print(f"\nArea stats before filter:\n{gdf['_area'].describe()}")
n_small = (gdf["_area"] < 1.0).sum()
gdf = gdf[gdf["_area"] >= 1.0].copy()
print(f"Removed {n_small} features with area < 1 m²")
# ---------------------------------------------------------------------------
# 5. Fix invalid geometries (buffer(0) on the remaining features)
# ---------------------------------------------------------------------------
n_invalid = (~gdf.geometry.is_valid).sum()
if n_invalid:
print(f"Fixing {n_invalid} invalid geometries with make_valid")
gdf.geometry = gdf.geometry.make_valid()
# After make_valid, some invalids might produce GeometryCollection etc.
# Keep only Polygon / MultiPolygon types
def ensure_poly(g):
if g.geom_type in ("Polygon", "MultiPolygon"):
return g
# GeometryCollection – extract Polygon/MultiPolygon parts
from shapely import GeometryCollection
if g.geom_type == "GeometryCollection":
parts = [p for p in g.geoms if p.geom_type in ("Polygon", "MultiPolygon")]
if not parts:
return None
from shapely.ops import unary_union
return unary_union(parts)
return None
orig = len(gdf)
gdf.geometry = gdf.geometry.apply(ensure_poly)
gdf = gdf[gdf.geometry.notna()].copy()
print(f"Dropped {orig - len(gdf)} features that had no valid poly after fix")
# ---------------------------------------------------------------------------
# 6. Convert all geometries to MultiPolygon
# (every feature must be MultiPolygon)
# ---------------------------------------------------------------------------
def to_multipoly(g):
if g.geom_type == "Polygon":
return MultiPolygon([g])
return g # already MultiPolygon
gdf.geometry = gdf.geometry.apply(to_multipoly)
assert gdf.geometry.geom_type.unique().tolist() == ["MultiPolygon"], \
"Not all geometries are MultiPolygon!"
# ---------------------------------------------------------------------------
# 7. Drop duplicate geometries — keep the row with the smallest record_seq
# ---------------------------------------------------------------------------
gdf["_wkt"] = gdf.geometry.apply(to_wkt)
# Sort by record_seq so the first occurrence per geometry is the earliest
gdf.sort_values("record_seq", inplace=True)
n_before_dedup = len(gdf)
gdf = gdf.drop_duplicates(subset="_wkt", keep="first").copy()
print(f"Dropped {n_before_dedup - len(gdf)} duplicate-geometry records (kept earliest record_seq)")
# ---------------------------------------------------------------------------
# 8. Recalculate area_m2 from the actual geometry
# ---------------------------------------------------------------------------
gdf["area_m2"] = gdf.geometry.area.round(2)
# ---------------------------------------------------------------------------
# 9. Select output columns
# ---------------------------------------------------------------------------
gdf_out = gdf[["parcel_id", "parcel_class", "district", "area_m2", "geometry"]].copy()
print(f"\nFinal record count: {len(gdf_out)}")
print(f"Columns: {list(gdf_out.columns)}")
print(f"Geometry types: {gdf_out.geometry.geom_type.value_counts().to_dict()}")
print(f"Any null geom: {gdf_out.geometry.isna().any()}")
print(f"Any empty geom: {gdf_out.geometry.is_empty.any()}")
print(f"All valid: {gdf_out.geometry.is_valid.all()}")
print(f"Area range: {gdf_out['area_m2'].min():.2f} – {gdf_out['area_m2'].max():.2f}")
# Check no duplicate geometries in output
dups = gdf_out.geometry.duplicated()
print(f"Duplicate geometries in output: {dups.sum()}")
# ---------------------------------------------------------------------------
# 10. Write to GeoParquet
# ---------------------------------------------------------------------------
gdf_out.to_parquet("/work/parcels_canonical.geoparquet", index=False)
print("\n✅ Wrote /work/parcels_canonical.geoparquet")Wrote /work/solve.py (5891 chars).
0:52
Let me run this script and see how it goes.
0:52
Bash
cd /work && python solve.py 2>&1
Total input features: 290
After null/empty removal: 290
CRS: EPSG:22992
Area stats before filter:
count 290.00000
mean 1001.43000
std 458.37842
min 0.00000
25% 1200.00000
50% 1200.00000
75% 1200.00000
max 1280.00000
Name: _area, dtype: float64
Removed 50 features with area < 1 m²
Dropped 50 duplicate-geometry records (kept earliest record_seq)
Final record count: 190
Columns: ['parcel_id', 'parcel_class', 'district', 'area_m2', 'geometry']
Geometry types: {'MultiPolygon': 190}
Any null geom: False
Any empty geom: False
All valid: True
Area range: 1200.00 – 1280.00
Duplicate geometries in output: 0
✅ Wrote /work/parcels_canonical.geoparquet
1:02
The output looks good. Let me verify the result by reading the geoparquet file to confirm everything is correct. Note: the count here is 290 features. 50 were removed for area < 1 m². And then 50 were dropped as duplicate geometries. Wait - but the 50 with area < 1m² are the same 20 invalid + 30 others with small area? Let me think... Actually, looking at the initial analysis: - 50 duplicate geometries (25 pairs) - 20 invalid geometries (all zero-area bowtie shapes) - 50 features with area < 1 m² So there's overlap. The 20 invalid (zero-area) ones are part of the 50 area < 1 m². And the 50 duplicate geometries dropped - those include the other members of each duplicate pair. Let me see: 290 total. 50 removed for area < 1m² = 240 remaining. Then 50 duplicates dropped from those 240 = 190 remaining. That seems right. But wait - we first remove small-area (<1 m²), which removes 50 features. The 50 area<1 features include the 20 invalid ones. Then from the remaining 240, duplicate geometries (by WKT) are removed - 50 duplicates found among the larger features. So 240 - 50 = 190 final. That seems correct. Let me verify the output file exists and sanity-check.
1:02
1:02
Bash
cd /work && python3 -c "
cd /work && python3 -c "
import geopandas as gpd
gdf = gpd.read_parquet('/work/parcels_canonical.geoparquet')
print('Shape:', gdf.shape)
print('Columns:', list(gdf.columns))
print('CRS:', gdf.crs)
print('dtypes:')
print(gdf.dtypes)
print()
print('Null geometries:', gdf.geometry.isna().sum())
print('Empty geometries:', gdf.geometry.is_empty.sum())
print('Invalid geometries:', (~gdf.geometry.is_valid).sum())
print('Non-MultiPolygon:', (gdf.geometry.geom_type != 'MultiPolygon').sum())
print('Duplicate geometries:', gdf.geometry.duplicated().sum())
print()
print('area_m2 range:', gdf['area_m2'].min(), '-', gdf['area_m2'].max())
print('area_m2 nulls:', gdf['area_m2'].isna().sum())
print()
print('parcel_id range:', gdf['parcel_id'].min(), '-', gdf['parcel_id'].max())
print('parcel_id nulls:', gdf['parcel_id'].isna().sum())
print('parcel_id duplicates:', gdf['parcel_id'].duplicated().sum())
print()
print('parcel_class distinct:', gdf['parcel_class'].unique())
print('district distinct:', gdf['district'].unique())
print()
# Check that all polygons are valid
print('Sample record:')
print(gdf.head(2).to_string())
" 2>&1
Shape: (190, 5)
Columns: ['parcel_id', 'parcel_class', 'district', 'area_m2', 'geometry']
CRS: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "ProjectedCRS", "name": "Egypt 1907 / Red Belt", "base_crs": {"name": "Egypt 1907", "datum": {"type": "GeodeticReferenceFrame", "name": "Egypt 1907", "ellipsoid": {"name": "Helmert 1906", "semi_major_axis": 6378200, "inverse_flattening": 298.3}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "id": {"authority": "EPSG", "code": 4229}}, "conversion": {"name": "Egypt Red Belt", "method": {"name": "Transverse Mercator", "id": {"authority": "EPSG", "code": 9807}}, "parameters": [{"name": "Latitude of natural origin", "value": 30, "unit": "degree", "id": {"authority": "EPSG", "code": 8801}}, {"name": "Longitude of natural origin", "value": 31, "unit": "degree", "id": {"authority": "EPSG", "code": 8802}}, {"name": "Scale factor at natural origin", "value": 1, "unit": "unity", "id": {"authority": "EPSG", "code": 8805}}, {"name": "False easting", "value": 615000, "unit": "metre", "id": {"authority": "EPSG", "code": 8806}}, {"name": "False northing", "value": 810000, "unit": "metre", "id": {"authority": "EPSG", "code": 8807}}]}, "coordinate_system": {"subtype": "Cartesian", "axis": [{"name": "Easting", "abbreviation": "E", "direction": "east", "unit": "metre"}, {"name": "Northing", "abbreviation": "N", "direction": "north", "unit": "metre"}]}, "scope": "Engineering survey, topographic mapping.", "area": "Egypt - onshore between 29\u00b0E and 33\u00b0E, offshore Mediterranean east of 29\u00b0E and offshore Gulf of Suez.", "bbox": {"south_latitude": 21.99, "west_longitude": 29, "north_latitude": 33.82, "east_longitude": 34.27}, "id": {"authority": "EPSG", "code": 22992}}
dtypes:
parcel_id int64
parcel_class str
district str
area_m2 float64
geometry geometry
dtype: object
Null geometries: 0
Empty geometries: 0
Invalid geometries: 0
Non-MultiPolygon: 0
Duplicate geometries: 0
area_m2 range: 1200.0 - 1280.0
area_m2 nulls: 0
parcel_id range: 1 - 210
parcel_id nulls: 0
parcel_id duplicates: 0
parcel_class distinct: <ArrowStringArray>
['commercial', 'industrial', 'agricultural', 'residential']
Length: 4, dtype: str
district distinct: <ArrowStringArray>
['Giza-East', 'Qalyubia-South', 'Cairo-Central']
Length: 3, dtype: str
Sample record:
parcel_id parcel_class district area_m2 geometry
0 1 commercial Giza-East 1200.0 MULTIPOLYGON (((640000 815000, 640030 815000, 640030 815040, 640000 815040, 640000 815000)))
1 2 industrial Qalyubia-South 1200.0 MULTIPOLYGON (((640050 815000, 640080 815000, 640080 815040, 640050 815040, 640050 815000)))
1:08
Everything looks good. Let me verify the file exists and also double-check the coordinate magnitudes are sensible for EPSG:22992 (meters). The coordinates ~640,000 / 815,000 are in meters, consistent with the Egypt Red Belt projection. The CRS is EPSG:22992 which is correct. Let me also check if there's any issues with the area_m2 precision - the original had area_m2 as 1200.0 which seems to match. And we have no nulls, no empties, no invalids, all MultiPolygon, no duplicate geometries, no duplicate parcel_ids. Let me confirm the output file exists.
1:08
1:08
Bash
ls -lh /work/parcels_canonical.geoparquet
-rw-r--r-- 1 agent agent 18K Jun 17 23:45 /work/parcels_canonical.geoparquet
1:15
Everything looks good. Let me summarize what was done.
1:15