7

I have a binary image with pixel resolution in units of degrees (EPSG:4326):

enter image description here

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:

enter image description here

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.]

Improve this question
3
  • 1
    This looks pretty useful: gist.github.com/habibutsu/8bbcc202a915e965c6a6d4f561d0e482
    – snowman2
    Sep 7, 2020 at 17:10
  • 1
    FYI, gdal_proximity can return distance in georeferenced units using the "-distunits GEO" option.
    – Brent Edwards
    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

1 Answer 1

Reset to default
1

I would follow these steps:

  • 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)
  • Create new raster from points' value (based-on source raster properties)

Since there will be millions of points for high resolution raster, this method will be time-consuming.

Improve this answer
1
  • I guess I knew how to "brute force" it like this. But as you mention, this would be very time consuming. It's too bad I can't just alter the distance equation in an off-the-shelf distance transform...OpenCV offers a cv.CV_DIST_USER that should allow a customized distance formula, but it's apparently only usable with the C API and I'm using Python.
    – Jon
    Sep 14, 2020 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.