3

I am having an issue that I have not been able to troubleshoot (see below for code to replicate). I am calculating the area in meters using ST_Area_Spheroid, but I am getting NaN as the result. I don't get NaN when using ST_Area, and I am also able to calculate the area in QGIS, so I am not sure where I am going wrong. Does anybody have any idea what could be causing this issue? I made sure to fix any invalid geometries using PostGIS's ST_MakeValid function.

import duckdb con = duckdb.connect() con.install_extension("spatial") con.load_extension("spatial") con.execute(""" DROP TABLE IF EXISTS geo_data; CREATE TABLE geo_data AS SELECT da_dguid, ST_Area_Spheroid(geom) as st_area_spheroid, ST_Area(geom) as st_area, geom FROM 'https://data-01.dataforcanada.org/processed/statistics_canada/boundaries/2021/digital_boundary_files/da_2021.parquet' WHERE isnan(ST_Area_Spheroid(geom)) """) issues = con.sql(""" SELECT * EXCLUDE geom FROM geo_data; """).df() issues
Improve this question
2
  • 1
    ST_Area_Spheroid assumes lat/long axis order. I have see it return NaN if the coordinates are passed in as long/lat. Try always_xy := true and see if that changes anything. Commented May 28 at 14:51
  • The only function that I see uses the always_xy parameter seems to be ST_Transform. Are you telling me to transform the geometry from EPSG:4326 to EPSG:4326? Commented May 28 at 18:23

1 Answer 1

Reset to default
4

There is a known, open issue with DuckDB's spatial extension when it comes to reading coordinates from files: Lat/long order when reading files · Issue #63 · duckdb/duckdb-spatial · GitHub. I encourage people to read through the issue discussion on GitHub.

Since the coordinates are already in EPSG:4326, the ST_FlipCoordinates function seems to handle the issue.

import duckdb con = duckdb.connect() con.install_extension("spatial") con.load_extension("spatial") con.execute(""" DROP TABLE IF EXISTS geo_data; CREATE TABLE geo_data AS SELECT da_dguid, ST_Area_Spheroid(st_flipcoordinates(geom)) as st_area_spheroid, ST_Area(st_flipcoordinates(geom)) as st_area, geom FROM 'https://data.dataforcanada.org/processed/statistics_canada/boundaries/2021/digital_boundary_files/da_2021.parquet' WHERE isnan(ST_Area_Spheroid(geom)) """) issues = con.sql(""" SELECT * EXCLUDE geom FROM geo_data; """).df() issues da_dguid st_area_spheroid st_area 0 2021S051259550180 1.006080e+05 0.000015 1 2021S051259550181 1.397088e+05 0.000020 2 2021S051259550178 1.085036e+05 0.000016 3 2021S051259550169 5.619264e+08 0.080906 4 2021S051259550182 1.631900e+06 0.000236 ... ... ... ... 19323 2021S051248140068 9.170059e+04 0.000012 19324 2021S051248142195 2.672939e+06 0.000361 19325 2021S051259530257 2.689198e+09 0.359978 19326 2021S051259530258 1.806210e+09 0.243301 19327 2021S051259530259 5.936013e+09 0.789738 [19328 rows x 3 columns]

I can't say whether the values are correct, but there are values instead of NaN.

Improve this answer
1
  • Yep, this did it! I confirmed the calculated values with QGIS and the values are identical. Commented May 30 at 15:09

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.