Search code examples
pythondjangogispostgisgeodjango

Can you do Azimuthal Equidistant projections natively in GeoDjango?

I am working on converting a small project I wrote to find overlapping boundaries of a shape file within a radius of a certain point. This original project was a mock up project I wrote using Shapely and GeoPandas, to make this more suitable for production, I am converting it all to GeoDjango.

There is one thing that is vital to this program, which is to create an equidistant projection of a circle on a map. I was able to do this with shapely objects using pyproj and functools.

Let it be known that this solution was found on stackoverflow and is not my original solution.

from shapely import geometry
from functools import partial

def createGeoCircle(lat, lng, mi):
    proj_wgs84 = pyproj.Proj(init='epsg:4326')
    aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lng} +x_0=0 +y_0=0'
    project = partial(
        pyproj.transform,
        pyproj.Proj(aeqd_proj.format(lat=lat, lng=lng)),
        proj_wgs84)
    buf = geometry.Point(0, 0).buffer(mi * 1.60934 * 1000)

    circle = transform(project, buf)
    return circle

I attempted to again use this solution and create a geoDjango MultiPolygon object from the shapely object, but it results in incorrect placement and shapes.

Here is the code I use to cast the shapely object coming from the above function.

shape_model(geometry=geos.MultiPolygon(geos.GEOSGeometry(createGeoCircle(41.378397, -81.2446768, 1).wkt)), state="CircleTest").save()

Here is the output in Django Admin. this picture is zoomed in to show the shape, but the location is in the middle of Antarctica. The coordinates given were meant to show in Ohio.

django admin geo output

To clear a few things up, my model is as follows:

class shape_model(geo_models.Model):
    state    = geo_models.CharField('State Territory ID', max_length=80)
    aFactor  = geo_models.FloatField()
    bFactor  = geo_models.FloatField()
    geometry = geo_models.MultiPolygonField(srid=4326)

I can get the location correct by simply using a geodjango point and buffer, but it shows up as an oval as it is not equidistant. If anyone has any suggestions or hints, I would be very appreciative to hear them!

Solution

  • Okay, I have found a solution to this problem. I used the shapely equidistant projection code and expanded it to convert it back to EPSG:4326. The updated function is as follows:

    def createGeoCircle(lat, lng, mi):
    point = geometry.Point(lat, lng)

    local_azimuthal_projection = f"+proj=aeqd +lat_0={lat} +lon_0={lng} +x_0=0 +y_0=0"
    proj_wgs84 = pyproj.Proj('epsg:4326')

    wgs84_to_aeqd = partial(
        pyproj.transform,
        proj_wgs84,
        pyproj.Proj(local_azimuthal_projection),
    )

    aeqd_to_wgs84 = partial(
        pyproj.transform,
        pyproj.Proj(local_azimuthal_projection),
        proj_wgs84,
    )

    point_transformed = transform(wgs84_to_aeqd, point)

    buffer = point_transformed.buffer(mi * 1.60934 * 1000)

    buffer_wgs84 = transform(aeqd_to_wgs84, buffer)

    return json.dumps(geometry.mapping(buffer_wgs84))

    I also dump the geometry mapping from this function so it can now be loaded directly into the geos MultiPolygon rather than using the wkt of the object. I load the circle into a model and save it using the following:

    shape_model(geometry=geos.MultiPolygon(geos.GEOSGeometry(createGeoCircle(41.378397, -81.2446768, 1))), state="CircleTest", aFactor=1.0, bFactor=1.0).save()

    FYI this is not a native geodjango solution and relies on many other packages. If someone has a native solution, I would greatly prefer that!