# Area interpolation between Shapefile (country borders) and WKT (multi)polygon elements in a CSV file
**URL:** https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458
**Category:** Geo
**Created:** [October 27, 2025, 1:26pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458 "2025-10-27T13:26:46Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [October 27, 2025, 1:26pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/1 "2025-10-27T13:26:46Z")
I have a [shapefile of administrative country borders](https://github.com/wmgeolab/geoBoundaries/raw/main/releaseData/CGAZ/geoBoundariesCGAZ_ADM0.zip) and a CSV file where one field is “geometry” and embed polygon and multipolygon data as "Well Known Text (WKT). Both are in EPSG:4326.
How can I retrieve for each element in the CSV file, the area under the various countries? I saw there is a package [AreaInterpolation.jl](https://github.com/JuliaGeo/AreaInterpolation.jl), but I don’t know how to get the data in the format required by that package.
This is the CSV file:
```julia-auto
Foo,Geometry
aaa,"MULTIPOLYGON (((10.992 11.493, 10.743 11.493, 10.743 11.743, 10.992 11.743, 10.992 11.493)), ((11.242 10.744, 10.992 10.744, 10.992 10.994, 10.743 10.994, 10.743 11.244, 10.992 11.244, 10.992 11.493, 11.242 11.493, 11.242 10.744)))"
bbb,"POLYGON ((11.161 10.799, 11.157 10.799, 11.157 10.803, 11.152 10.803, 11.152 10.808, 11.161 10.808, 11.161 10.812, 11.161 10.821, 11.166 10.821, 11.174 10.821, 11.179 10.821, 11.183 10.821, 11.183 10.816, 11.188 10.816, 11.188 10.812, 11.183 10.812, 11.183 10.808, 11.183 10.803, 11.188 10.803, 11.188 10.794, 11.183 10.794, 11.183 10.79, 11.179 10.79, 11.174 10.79, 11.174 10.794, 11.17 10.794, 11.17 10.799, 11.166 10.799, 11.161 10.799))"
ccc,"POLYGON ((4.497 13.49, 4.497 12.991, 3.997 12.991, 3.997 13.24, 3.498 13.24, 3.498 13.49, 4.497 13.49))"
```
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [October 27, 2025, 2:21pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/2 "2025-10-27T14:21:12Z")
```juila
using WellKnownGeometry, GeoFormatTypes, GeometryOps, DataFrames, Proj
GeometryOps.area.((Geodesic(),), WellKnownText.(DataFrame(csv).Geometry))
```
---
### Author:  [@evetion](https://discourse.julialang.org/u/evetion)
#### Post date: [October 27, 2025, 2:21pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/3 "2025-10-27T14:21:15Z")
I’m not sure how you want your output to look like, as you want the area for each country per your csv polygons, but here’s an example of the total area for each polygon.
This uses GDAL underneath, but we reproject the geometries first to an equal area projection ([WGS 84 / NSIDC EASE-Grid 2.0 Global - EPSG:6933](https://epsg.io/6933)), so the area will be in metres squared. Note that this is a bit inefficient, as we do an intersects on all countries three times. For a larger csv you’d want to setup a spatial index first.
```julia
using GeoDataFrames
import GeoInterface as GI
countries = GeoDataFrames.read("geoBoundariesCGAZ_ADM0.shp")
countries = reproject(countries, EPSG(4326), EPSG(6933))
# rename Geometry column to wkt in the csv file for automatic detection
# otherwise we need to pass GDAL options specifically
polygons = GeoDataFrames.read("test.csv")
polygons = reproject(polygons, EPSG(4326), EPSG(6933))
polygons.countries_area_km2 = map(polygons.geometry) do geom
# Find intersecting countries (faster than calculating intersection)
icountries = countries.geometry[GI.intersects.(countries.geometry, (geom,))]
# For intersecting countries, calculate the actual intersection
intersections = GI.intersection.(icountries, (geom,))
sum(GI.area.(intersections) / 1e6) # in km²
end
3-element Vector{Float64}:
3765.9015712420355
8.600586983524826
4486.604217175244
```
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [October 27, 2025, 2:23pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/4 "2025-10-27T14:23:23Z")
Isn’t using `GeometryOps.area` much easier than that?
I though GI.area was deprecated
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [October 27, 2025, 2:24pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/5 "2025-10-27T14:24:35Z")
Why do you reproject both data and administrative borders to EPSG 6933 ? Because otherwise there is a bias in the original EPSG and the computation of the are is not correct ?
---
### Author:  [@evetion](https://discourse.julialang.org/u/evetion)
#### Post date: [October 27, 2025, 2:26pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/6 "2025-10-27T14:26:25Z")
Your input data is in degrees (curved polygons on the sphere), while most computations assume planar geometries. So not only would the intersections be slightly wrong, your area results would be in square degrees, which is probably not what you want.
---
### Author:  [@evetion](https://discourse.julialang.org/u/evetion)
#### Post date: [October 27, 2025, 2:27pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/7 "2025-10-27T14:27:30Z")
Yes, ideally we use GeometryOps here, but it has no geodetic area yet, and I found a bug with intersecting geometries.
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [October 27, 2025, 2:28pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/8 "2025-10-27T14:28:10Z")
It does have geodesic area, just need to load Proj
---
### Author:  [@evetion](https://discourse.julialang.org/u/evetion)
#### Post date: [October 27, 2025, 2:44pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/9 "2025-10-27T14:44:04Z")
I stand corrected. That said, while `area` has a `Geodesic` trait, `intersects`/`intersection` doesn’t AFAIK, so it can be quite wrong (although probably acceptable for the small scale, low latitude example here). Anyway, using GeometryOps is mostly replacing GeoInterface with GeometryOps. Here’s the code, but it fails on the intersections, maybe I’m doing something wrong. I’ll make an issue @asinghvi17
```julia
using CSV
using WellKnownGeometry
using GeoFormatTypes
import GeometryOps as GO
import Proj
using GeoDataFrames
import GeoInterface as GI
countries = GeoDataFrames.read(
"geoBoundariesCGAZ_ADM0.shp",
)
polygons = CSV.read("test.csv", DataFrame)
polygons.geometry =
GeoFormatTypes.WellKnownText.((GeoFormatTypes.Geom(),), polygons.Geometry)
polygons.countries_area_km2 = map(polygons.geometry) do geom
icountries = countries.geometry[GO.intersects.(countries.geometry, (geom,))]
intersections = GO.intersection.(icountries, (geom,)) # fails here, also with setting with `target` kwarg
sum(GO.area.((Geodesic(),), intersections) / 1e6) # in km²
end
```
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [October 27, 2025, 3:28pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/10 "2025-10-27T15:28:53Z")
I have some issues with this script (different package versions ??)
```julia-auto
countries = GeoDataFrames.read("adm_boundaries/geoBoundariesCGAZ_ADM0.shp")
countries = reproject(countries, EPSG(4326), EPSG(6933))
```
This returns me an error `ERROR: GDALError (CE_Failure, code 1): PROJ: cea: Invalid latitude`
Strange enough, it works if I apply it in a loop:
```julia-auto
countries = GeoDataFrames.read("adm_boundaries/geoBoundariesCGAZ_ADM0.shp")
countries.geometry = [reproject(countries[i,"geometry"], EPSG(4326), EPSG(6933)) for i in 1:size(countries,1)]
```
For the data file, I have the issue that `reproject` leaves the Geometry field as a (unchanged) text field, likely related to the fact that I have the wkt field with a different name. I solved similarly with a loop and using ArcGDAL.
```julia-auto
import ArchGDAL
polygons = GeoDataFrames.read("test.csv")
polygons.geometry = [reproject(ArchGDAL.fromWKT(polygons[i,"Geometry"]), EPSG(4326), EPSG(6933)) for i in 1:size(polygons,1)]
```
So now I have both countries and polygons geometries as vector of reprojected polygons `Vector{GeoInterface.Wrappers.WrapperGeometry{false, false, T, EPSG{1}} where T}`
Still, the `intersects` call fail with error:
```julia-auto
julia> polygons.countries_area_km2 = map(polygons.geometry) do geom
# Find intersecting countries (faster than calculating intersection)
icountries = countries.geometry[GI.intersects.(countries.geometry, (geom,))]
# For intersecting countries, calculate the actual intersection
intersections = GI.intersection.(icountries, (geom,))
sum(GI.area.(intersections) / 1e6) # in km²
end
ERROR: MethodError: no method matching intersects(::PolygonTrait, ::MultiPolygonTrait, ::GeoInterface.Wrappers.Polygon{…}, ::GeoInterface.Wrappers.MultiPolygon{…})
The function `intersects` exists, but no method is defined for this combination of argument types.
Closest candidates are:
intersects(::Any, ::Any)
@ GeoInterface ~/.julia/packages/GeoInterface/X1jn3/src/interface.jl:501
intersects(::Union{CircularStringTrait, CompoundCurveTrait, CurvePolygonTrait, GeometryCollectionTrait, LineStringTrait, LinearRingTrait, MultiLineStringTrait, MultiPointTrait, MultiPolygonTrait, MultiSurfaceTrait, PointTrait, PolygonTrait, PolyhedralSurfaceTrait, TINTrait, TriangleTrait}, ::Union{CircularStringTrait, CompoundCurveTrait, CurvePolygonTrait, GeometryCollectionTrait, LineStringTrait, LinearRingTrait, MultiLineStringTrait, MultiPointTrait, MultiPolygonTrait, MultiSurfaceTrait, PointTrait, PolygonTrait, PolyhedralSurfaceTrait, TINTrait, TriangleTrait}, ::ArchGDAL.AbstractGeometry, ::ArchGDAL.AbstractGeometry)
@ ArchGDAL ~/.julia/packages/ArchGDAL/N2rQ4/src/geointerface.jl:198
Stacktrace:
[1] intersects(a::GeoInterface.Wrappers.Polygon{…}, b::GeoInterface.Wrappers.MultiPolygon{…})
@ GeoInterface ~/.julia/packages/GeoInterface/X1jn3/src/interface.jl:501
[2] _broadcast_getindex_evalf
@ ./broadcast.jl:678 [inlined]
[3] _broadcast_getindex
@ ./broadcast.jl:651 [inlined]
[4] getindex
@ ./broadcast.jl:610 [inlined]
[5] copy
@ ./broadcast.jl:911 [inlined]
[6] materialize(bc::Base.Broadcast.Broadcasted{…})
@ Base.Broadcast ./broadcast.jl:872
[7] (::var"#23#24")(geom::GeoInterface.Wrappers.MultiPolygon{false, false, Vector{…}, Nothing, EPSG{…}})
@ Main ./REPL[6]:3
[8] iterate
@ ./generator.jl:48 [inlined]
[9] _collect(c::Vector{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
@ Base ./array.jl:811
[10] collect_similar(cont::Vector{…}, itr::Base.Generator{…})
@ Base ./array.jl:720
[11] map(f::Function, A::Vector{GeoInterface.Wrappers.WrapperGeometry{false, false, T, EPSG{1}} where T})
@ Base ./abstractarray.jl:3371
[12] top-level scope
@ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.
```
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [October 27, 2025, 3:46pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/11 "2025-10-27T15:46:13Z")
I can get `intersect`/`intersection` working by postponing the reprojection after the intersection (shoudn’t be an issue, right?), but I still have issue computing the area on the reprojected intersection.
My intersection is now a GeoInterface _wrapper_:
```julia-auto
julia> c1_d1_intersection_reprojection = GeoDataFrames.reproject(c1_d1_intersection, GeoDataFrames.EPSG(4326), GeoDataFrames.EPSG(6933))
GeoInterface.Wrappers.MultiPolygon{false, false}([GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.0365521087353812e6, 1.4566451802948604e6), … (3) … , (1.0365521087353812e6, 1.4566451802948604e6)], crs = "EPSG:6933")], crs = "EPSG:6933"), GeoInterface.Wrappers.Polygon([GeoInterface.Wrappers.LinearRing([(1.0605771925178545e6, 1.3628402165805716e6), … (7) … , (1.0605771925178545e6, 1.3628402165805716e6)], crs = "EPSG:6933")], crs = "EPSG:6933")], crs = "EPSG:6933")
```
How do I convert this wrapper of a (multi)polygon to an actual (multi)polygon with the “reprojected” coordinates ?
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [October 27, 2025, 4:14pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/12 "2025-10-27T16:14:29Z")
SOLVED using `ArcGDAL.reproject` instead of `GeoDataFrames.reproject` (that creates a wrapper).
Reprojecting before or after the intersection provide the same result:
```julia-auto
import Proj
import GeoDataFrames
import GeoFormatTypes
import ArchGDAL
import GeoInterface
map_path = joinpath(@ __DIR__ ,"adm_boundaries","geoBoundariesCGAZ_ADM0.shp")
data_path = "testdata.csv"
map = GeoDataFrames.read(map_path)
map = map[map.shapeType .== "ADM0", :] # filter only countries
data = GeoDataFrames.read(joinpath(@ __DIR__ ,"testdata.csv"))
# Testing with one geometry
d1 = ArchGDAL.fromWKT(data[1,"Geometry"])
c1 = map[map.shapeGroup .== "NGA","geometry"][1]
c1_d1_intersection = GeoInterface.intersection(c1,d1)
c1_d1_intersection_reprojection = ArchGDAL.reproject(c1_d1_intersection, GeoFormatTypes.EPSG(4326), GeoFormatTypes.EPSG(6933))
a1 = GeoInterface.area(c1_d1_intersection_reprojection) / 1e6
d1_rep = ArchGDAL.reproject(d1, GeoFormatTypes.EPSG(4326), GeoFormatTypes.EPSG(6933))
c1_rep = ArchGDAL.reproject(c1, GeoFormatTypes.EPSG(4326), GeoFormatTypes.EPSG(6933))
c1_d1_intersection2 = ArchGDAL.intersection(c1,d1)
a2 = GeoInterface.area(c1_d1_intersection2) / 1e6
@assert a1 ≈ a2
```
---
### Author:  [@evetion](https://discourse.julialang.org/u/evetion)
#### Post date: [October 27, 2025, 8:21pm UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/13 "2025-10-27T20:21:26Z")
I’m happy to hear you solved your problem 🙂, and good to see there’s no real difference (in this case) between the two areas when reprojecting before or after.
It can be confusing to use different level of abstractions as you’re doing here. And you’re right about different package versions, I might’ve been on the main branch of GeoDataFrames when I first tried it. I’m not sure about your manual loops, it might be using different methods (given individual elements versus the vector of elements in a DataFrame). Some pointers:
Reading the CSV with the `Geometry` column (not renamed to `wkt` as I proposed) can be done with GeoDataFrames like so (using [Comma Separated Value (.csv) — GDAL documentation](https://gdal.org/en/stable/drivers/vector/csv.html) as reference), otherwise the Geometry is not parsed (but stays a String as you noticed).
```julia
polygons = GeoDataFrames.read("test.csv", options=["GEOM_POSSIBLE_NAMES=Geometry", "KEEP_GEOM_COLUMNS=NO"])
```
The invalid latitude error is fixed on the main branch, soon to be released, currently you want to use
```julia
countries = reproject(countries, EPSG(4326), EPSG(6933); always_xy=true) # forces longitude, latitude order
```
Finally, there’s indeed two ways of doing this, either with (Arch)GDAL (using non-wrapped geometries everywhere though), or using GeometryOps (but I hit an error doing the latter). Mixing them doesn’t really work.
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [October 28, 2025, 3:11am UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/14 "2025-10-28T03:11:52Z")
I think we need to stop recommending people use any of these `GeoInterface.intersection` etc methods, they’re so fragile because only ArchGDAL implemented them.
`GeometryOps.intersection` (etc) works on anything and is faster.
`GeometryOps.reproject`is also better than `ArchGDAL.reproject` (and I wrote `ArchGDAL.reproject`)
Its unfortunate that there was a bug in GeometryOps because it really would take a fraction of this code to do the same thing.
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [October 28, 2025, 3:14am UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/15 "2025-10-28T03:14:31Z")
The key point in “replacing GeoInterface with GeometryOps” is that the GeometryOps code will work on literally any geometry and the GeoInterface code wont.
---
### Author:  [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [October 30, 2025, 1:01am UTC](https://discourse.julialang.org/t/area-interpolation-between-shapefile-country-borders-and-wkt-multi-polygon-elements-in-a-csv-file/133458/16 "2025-10-30T01:01:00Z")
> [@sylvaticus](#):
>
> `julia> c1_d1_intersection_reprojection = GeoDataFrames.reproject(c1`
The coordinates look to be in the CRS you want them to be in, though? The “Wrapper” thing can be a bit misleading but this is a real, instantiated geometry with those coordinates, it’s not lazy.