Create Polygon with Length, Width, and a center coordinate

Question

Using Shapely or like library, How can I create a rectangular polygon with known Length, Width (in meters)(Drawing over EPSG: 3857), and center coordinates (EPSG: 4326)? The output needs to be in GeoJSON format.

import shapely.geometry
from shapely import affinity
from shapely.geometry import shape, LineString
import geopandas as gp
import geojson
import json
from ast import literal_eval
from functools import partial
import pyproj
from shapely.ops import transform


def main():
    # create geodataframe with some variables
    gf = gp.GeoDataFrame({'lat':38.6159031816694, 'lon':-90.2106388, 'width':100, 'height':60}, index=[0])

    # set default crs for dataframe
    gf.crs = {'init' :'epsg:4326'}

    # create center as a shapely geometry point type and set geometry of dataframe to this
    gf['center'] = gf.apply(lambda x: shapely.geometry.Point(x['lon'], x['lat']), axis=1)
    gf = gf.set_geometry('center')

    # change crs of dataframe to projected crs to enable use of distance for width/height
    gf = gf.to_crs(epsg = 3857)

    # create polygon using width and height
    gf['center'] = gf['center'].buffer(1)
    gf['polygon'] = gf.apply(lambda x: shapely.affinity.scale(x['center'], x['width'], x['height']), axis=1)
    gf = gf.set_geometry('polygon')
    # gf = affinity.rotate(gf, 11.60, origin='centroid')

    # return to inital input crs
    # gf = gf.to_crs(epsg = 4326)

    geopoly = gf['polygon'].to_json()
    a = geojson.loads(geopoly)
    # print(a.bbox)
    startbbox(a.bbox)


def bbox2wkt(minx, miny, maxx, maxy):
# first validate bbox values
    assert isinstance(minx,float) or isinstance(minx,int)
    assert isinstance(miny,float) or isinstance(miny,int)
    assert isinstance(maxx,float) or isinstance(maxx,int)
    assert isinstance(maxy,float) or isinstance(maxy,int)
    assert (minx < maxx), 'failed: %s is not < %s' % (minx,maxx)
    assert (miny < maxy), 'failed: %s is not < %s' % (miny,maxy)
    return "[%(minx)s, %(miny)s], [%(minx)s, %(maxy)s], [%(maxx)s, %(maxy)s], " \
           "[%(maxx)s, %(miny)s], [%(minx)s, %(miny)s]" % locals()


def startbbox(str1):
    # bbox = '-90.21153711528392,38.61548204271258,-90.20974048471568,38.6163243181534'
    bbox = map(str, str1)
    bbox = ','.join(bbox)
    rtn = bbox2wkt(*map(float, bbox.split(',')))
    data = list(literal_eval(rtn))
    gd_info = dict(type="Polygon", coordinates=[data])
    s = json.dumps(gd_info)
    g1 = geojson.loads(s)
    g2 = shape(g1)
    gf = affinity.rotate(g2, 90, origin='centroid')
    project = partial(
        pyproj.transform,
        pyproj.Proj(init='epsg:3857'),  # source coordinate system
        pyproj.Proj(init='epsg:4326'))  # destination coordinate system
    g2 = transform(project, gf)
    wow3 = geojson.dumps(g2)
    wow4 = json.loads(wow3)
    with open('data222.json', 'w') as f:
        geojson.dump(wow4, f)


if __name__ == "__main__":
    main()

Now I just need to clean it up.

Answer 1

Very much depends on what shape you're looking for. Assuming some sort of spherical object and not a rectangle or square:

import shapely
import geopandas as gp

#create geodataframe with some variables
gf = gp.GeoDataFrame({'lat':0, 'lon':0, 'width':100, 'height':20}, index=[0])

#set default crs for dataframe
gf.crs = {'init' :'epsg:4326'}

#create center as a shapely geometry point type and set geometry of dataframe to this
gf['center'] = gf.apply(lambda x: shapely.geometry.Point(x['lon'], x['lat']), axis=1)
gf = gf.set_geometry('center')

#change crs of dataframe to projected crs to enable use of distance for width/height
gf = gf.to_crs(epsg = 3857)


#create polygon using width and height
gf['center'] = gf['center'].buffer(1)
gf['polygon'] = gf.apply(lambda x: shapely.affinity.scale(x['center'], x['width'], x['height']), axis=1)
gf = gf.set_geometry('polygon')

#return to inital input crs
gf = gf.to_crs(epsg = 4326)

#plot the polygon
gf.plot()

Answer 2

@saph_top this helps:

def image_poly(imgar):
    polys = []
    for cent in imgar:
        lat = cent['coords'][1]
        lng = cent['coords'][0]
        prps = cent['props']
        wid= prps['Image_Width']
        hite = prps['Image_Height']
        head = prps['Heading']
        img_n = prps['File_Name']
        focal_lgth = prps['Focal_Length']
        alt = float(prps["Relative_Altitude"])
        calc = get_area(wid, hite, alt, focal_lgth)
        # polyg = polygen(lng, lat, calc[1], calc[0])

        # create geodataframe with some variables
        gf = gp.GeoDataFrame({'lat': lat, 'lon': lng, 'width': calc[0], 'height': calc[1]}, index=[1])

        gf.crs = fiona.crs.from_epsg(4326)

        # create center as a shapely geometry point type and set geometry of dataframe to this
        gf['center'] = gf.apply(lambda x: shapely.geometry.Point(x['lon'], x['lat']), axis=1)
        gf = gf.set_geometry('center')

        # change crs of dataframe to projected crs to enable use of distance for width/height
        gf = gf.to_crs(epsg = 3857)

        # create polygon using width and height
        gf['center'] = shapely.geometry.box(*gf['center'].buffer(1).total_bounds)
        gf['polygon'] = gf.apply(lambda x: shapely.affinity.scale(x['center'], x['width'], x['height']), axis=1)
        gf = gf.set_geometry('polygon')
        geopoly = gf['polygon'].to_json()
        g1 = geojson.loads(geopoly)
        gh = g1[0].geometry
        g2 = shape(gh)

        # Change Heading / orientation
        head_ck = 0
        if head > head_ck:
            header = 360 - float(head)
        else:
            header = float(head) - 360
        ngf = affinity.rotate(g2, header, origin='centroid')
        # Transform back to 4326
        project = partial(
            pyproj.transform,
            pyproj.Proj(init='epsg:3857'),  # source coordinate system
            pyproj.Proj(init='epsg:4326'))  # destination coordinate system
        g2 = transform(project, ngf)
        # Create GeoJSON
        wow3 = geojson.dumps(g2)
        wow4 = json.loads(wow3)
        propy = dict(image=img_n)
        gd_feat = dict(Feature="Polygon", properties=propy, geometry=wow4)
        polys.append(gd_feat)
    return polys


def get_area(wd, ht, alt, fl):
    sw = 9
    gsd = (sw * alt * 100) / (fl * wd)
    grd_wdth = ((gsd * wd) / 100)
    grd_ht = ((gsd * ht) / 100)
    return [grd_wdth, grd_ht]