Dutch Portolan Demo

# BRP Gewaspercelen (Agricultural Crop Parcels) — RVO / Netherlands

## What This Dataset Is

Every agricultural parcel in the Netherlands — 2,331,498 polygons — with its registered
crop type, from the Basisregistratie Gewaspercelen (BRP), the Dutch national crop parcel
registration. This is the 2025 definitive edition ("definitief"), finalized after
verification by RVO.

The BRP is an annual snapshot reflecting what farmers register each May 15, the deadline
for Common Agricultural Policy (CAP) subsidy applications. It is the authoritative record
of agricultural land use in the Netherlands: which crop is grown on which parcel, in which
year.

**Key concepts:**
- **category** = broad crop group (e.g., "Grasland", "Aardappelen", "Granen")
- **gewas** = specific crop within that group (e.g., "Blijvend grasland",
  "Consumptieaardappelen", "Wintertarwe")
- **gewascode** = numeric code for the specific crop (links to RVO's official crop code list)
- **definitief** = the edition is finalized after verification (vs. "voorlopig" = preliminary)

**Source:** https://service.pdok.nl/rvo/gewaspercelen/atom/downloads/brpgewaspercelen_definitief_2025.gpkg
**Provider:** RVO (Rijksdienst voor Ondernemend Nederland / Netherlands Enterprise Agency)
**License:** CC0 (public domain)

## How to Access

This is a GeoParquet file (1.7 GB) with geometries in EPSG:28992 (Amersfoort / RD New),
the Dutch national coordinate system. Use DuckDB with the spatial extension.

```python
import duckdb

con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;")

URL = 'https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet'

df = con.execute(f"""
    SELECT * FROM read_parquet('{URL}')
    LIMIT 5
""").df()
```

The file is 1.7 GB. DuckDB streams it efficiently using HTTP range requests — you don't
need to download the full file. For spatial queries, DuckDB can skip irrelevant row groups
using the Parquet metadata.

## Schema — Field Meanings

| Field | Type | Meaning |
|-------|------|---------|
| `geom` | WKB Polygon | Parcel boundary in **EPSG:28992** (Amersfoort / RD New). This is the Dutch national CRS — coordinates are in metres. |
| `id` | int64 | Feature ID. |
| `category` | string | **Broad crop category** — groups similar crops together (e.g., "Grasland" covers all grass types, "Aardappelen" covers all potato varieties). There are approximately 30-40 distinct categories. |
| `gewas` | string | **Specific crop name** within the category. More granular than category (e.g., under "Aardappelen": "Consumptieaardappelen", "Pootaardappelen", "Zetmeelaardappelen"). There are 300+ distinct gewas values. |
| `gewascode` | int32 | **Numeric crop code** corresponding to the gewas. Links to RVO's official crop code table. |
| `jaar` | int32 | Registration year — 2025 for this edition. |
| `status` | string | Always "definitief" in this edition (finalized after verification). |

## Category vs Gewas

The **category** field provides a high-level grouping useful for thematic mapping and
summary statistics. The **gewas** field gives the precise crop registration. Examples:

| category | gewas (examples) | gewascode |
|----------|-------------------|-----------|
| Grasland | Blijvend grasland, Tijdelijk grasland | 265, 266 |
| Aardappelen | Consumptieaardappelen, Pootaardappelen, Zetmeelaardappelen | 2014, 2015, 2016 |
| Granen | Wintertarwe, Zomertarwe, Wintergerst | 233, 234, 235 |
| Mais | Snijmais, Korrelsnijmais | 259, 260 |

Use `category` for broad analysis ("how much grassland vs arable?") and `gewas` for
specific crop statistics.

## Geometry Notes

- CRS is **EPSG:28992** (Amersfoort / RD New) — the standard Dutch national coordinate
  system. Coordinates are in metres, not degrees.
- To convert to WGS84 (lat/lon) for mapping, use `ST_Transform(ST_GeomFromWKB(geom), 'EPSG:28992', 'EPSG:4326')`.
- Geometry column is named `geom` (WKB encoded)
- All geometries are Polygons
- Bounding box in WGS84: approximately 3.21 to 7.24 longitude, 50.73 to 53.58 latitude
  (covers all agricultural land in the Netherlands)

## Useful Query Patterns

### Preview the data

```sql
SELECT id, category, gewas, gewascode, jaar, status
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet')
LIMIT 10
```

### Crop distribution — number of parcels per category

```sql
SELECT category, COUNT(*) AS parcels, ROUND(COUNT(*) * 100.0 / SUM(COUNT(*)) OVER (), 1) AS pct
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet')
GROUP BY category
ORDER BY parcels DESC
```

### Total area per crop category (in hectares)

```sql
INSTALL spatial; LOAD spatial;

SELECT category,
       COUNT(*) AS parcels,
       ROUND(SUM(ST_Area(ST_GeomFromWKB(geom))) / 10000, 0) AS area_ha
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet')
GROUP BY category
ORDER BY area_ha DESC
```

Note: since the CRS is EPSG:28992 (metres), `ST_Area` returns square metres directly — no
spheroidal calculation needed. Divide by 10,000 for hectares.

### Top 20 most common specific crops (gewas)

```sql
SELECT gewas, gewascode, category, COUNT(*) AS parcels
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet')
GROUP BY gewas, gewascode, category
ORDER BY parcels DESC
LIMIT 20
```

### Find parcels near a location (e.g., near Wageningen — 51.97, 5.67)

Since the data is in EPSG:28992, you need to either transform the point or use bbox columns
if available. Here is the approach using coordinate transformation:

```sql
INSTALL spatial; LOAD spatial;

-- Wageningen in RD New coordinates (approximately x=173600, y=443600)
SELECT id, category, gewas,
       ST_AsText(ST_GeomFromWKB(geom)) AS wkt
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet')
WHERE ST_DWithin(
    ST_GeomFromWKB(geom),
    ST_Point(173600, 443600),
    1000
)
```

This finds all parcels within 1 km of Wageningen. Because the CRS is in metres, the
distance parameter is in metres.

### Parcel size distribution

```sql
INSTALL spatial; LOAD spatial;

SELECT
    CASE
        WHEN ST_Area(ST_GeomFromWKB(geom)) / 10000 < 1 THEN '< 1 ha'
        WHEN ST_Area(ST_GeomFromWKB(geom)) / 10000 < 5 THEN '1-5 ha'
        WHEN ST_Area(ST_GeomFromWKB(geom)) / 10000 < 10 THEN '5-10 ha'
        WHEN ST_Area(ST_GeomFromWKB(geom)) / 10000 < 50 THEN '10-50 ha'
        ELSE '50+ ha'
    END AS size_class,
    COUNT(*) AS parcels
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet')
GROUP BY 1
ORDER BY 1
```

### Spatial join with municipalities (bestuurlijke gebieden)

To find crop distribution per municipality, join with the bestuurlijke_gebieden dataset:

```sql
INSTALL spatial; LOAD spatial;

SELECT g.naam AS municipality, b.category, COUNT(*) AS parcels
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/cbs/bestuurlijke_gebieden/gemeentegebied.parquet') g
JOIN read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet') b
  ON ST_Intersects(ST_GeomFromWKB(g.geom), ST_Centroid(ST_GeomFromWKB(b.geom)))
GROUP BY g.naam, b.category
ORDER BY g.naam, parcels DESC
```

Note: this is a heavy query on 2.3M parcels. Consider filtering to a specific municipality
or category first.

### Proximity to Natura 2000 areas

Find crop parcels within 500m of a Natura 2000 protected area:

```sql
INSTALL spatial; LOAD spatial;

SELECT b.id, b.category, b.gewas, n.naam_n2k AS natura2000_area
FROM read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/natura2000/natura2000.parquet') n
JOIN read_parquet('https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet') b
  ON ST_DWithin(ST_GeomFromWKB(n.geometry), ST_GeomFromWKB(b.geom), 500)
LIMIT 100
```

Both datasets are in EPSG:28992, so the 500m buffer works directly.

### Load into GeoPandas

```python
import geopandas as gpd

gdf = gpd.read_parquet(
    'https://data.source.coop/cholmes/portolan-nl/rvo/brp_gewaspercelen/brp_gewaspercelen.parquet',
    columns=['id', 'category', 'gewas', 'gewascode', 'geom']
)
# CRS will be EPSG:28992 — to convert to WGS84:
# gdf = gdf.to_crs(epsg=4326)
```

## Related Datasets in This Catalog

- **Natura 2000 Protected Areas** (`rvo/natura2000/`) — 162 protected nature areas. Useful
  for proximity analysis: which crops are grown near protected habitats?
- **Bestuurlijke Gebieden** (`cbs/bestuurlijke_gebieden/`) — municipality, province, and
  country boundaries. Spatial join to aggregate crop data by administrative unit.
- **National Parks** (`rvo/nationale_parken/`) — national park boundaries for similar
  proximity analysis.

## Use Cases

- **Agricultural land-use mapping:** What crops are grown where across the Netherlands?
- **CAP subsidy analysis:** The BRP underpins the EU Common Agricultural Policy payments —
  parcels must be registered to receive subsidies.
- **Environmental monitoring:** Track grassland vs arable ratios, crop diversity, and
  agricultural intensification trends over time (compare across years).
- **Nature impact assessment:** Which crops are grown near Natura 2000 areas? Nitrogen
  deposition from agriculture is a major concern for Dutch nature policy.
- **Food security:** Total hectares per crop type across the country.

## Caveats

- **EPSG:28992 coordinates:** Unlike some other datasets in this catalog, coordinates are
  in the Dutch RD New system (metres), not WGS84 (degrees). You must transform for web
  mapping or to join with WGS84 datasets.
- **Annual snapshot:** This is the state as of May 15, 2025. Crop registration may differ
  from what is actually grown (though discrepancies can lead to subsidy sanctions).
- **Large file:** 1.7 GB — DuckDB handles this efficiently with HTTP range requests, but
  downloading the full file will take time on slower connections.
- **No yield or economic data:** The BRP only records which crop is registered on each
  parcel, not yields, revenues, or farming practices.
- **Geometry column name:** The geometry column is named `geom`, not `geometry`. Adjust
  queries accordingly.

## Also Available As

- **PMTiles (vector tiles):** `brp_gewaspercelen.pmtiles` — for web map visualization with
  MapLibre GL JS or similar. Features are simplified/dropped at lower zoom levels.
- **Source GeoPackage:** Download from PDOK at
  https://service.pdok.nl/rvo/gewaspercelen/atom/downloads/brpgewaspercelen_definitief_2025.gpkg