# Haversine based distance transform

I have a binary image with pixel resolution in units of degrees (EPSG:4326): In this case, all the light-blue pixels are water (value = 1) and purple pixels are not water (value = 0). I would like to compute a distance transform of this image, where the result is each pixel's distance away from the nearest "on" (water) pixel. This is doable with scipy: However, the returned distances are Euclidean with respect to the row, column coordinates of each pixel. Does anyone know of a package or function that will compute a distance transform using the Haversine formula on the lon, lat coordinates rather than the row, col coordinates?

[I know I can reproject to a length-preserving CRS and multiply by the resolution.]

• This looks pretty useful: gist.github.com/habibutsu/8bbcc202a915e965c6a6d4f561d0e482 Sep 7, 2020 at 17:10
• FYI, gdal_proximity can return distance in georeferenced units using the "-distunits GEO" option. Sep 8, 2020 at 15:10
• @BrentEdwards Nice, I'd never seen that function before! I wonder if it really solves the problem though. E.g. I can return georeferenced units for 4326 (degrees), but the problem is still that a degree of latitude and degree of longitude don't represent the same distance, so I don't know how valid converting the output to distances using a scaling factor would be. I could apply a scaling factor based on latitude, but even that's not quite correct because the return of gdal_proximity lumps lat/lon degrees together in the output. Anyway, thanks for the heads up!
– Jon
Sep 8, 2020 at 20:58

• Create points from individual pixel's center, assign each pixel value and coordinate of its center to the corresponding point. Review this post. (`rasterio`, `geopandas`)
• Collect all water points to one multipoint object. (`geopandas`)
• Calculate haversine distance between a point and the multipoint and assign the distance to the point. In this step, the result is each point's distance away from the nearest point in the multipoint (water points). Review this post. (`geopandas`/`shapely`)