1

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.

  • Please Edit the question to include your best attempt to solve this programming question – Vince Feb 2 at 22:35
3

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()
  • Actually I am looking for a rectangle. My bad, I should have been more specific. – Dean Feb 3 at 0:44
  • 2
    @Dean, this saph_top's answer is too good just because you have not shown any attempt to solve your question. In addition to the omission that you want a rectangle polygon, knowing that the earth's surface is curved, you should clarify about over which plane coordinate system you want to draw it. – Gabriel De Luca Feb 3 at 1:34
  • @GabrielDeLuca - that additional information has been added. – Dean Feb 3 at 15:19
  • @Dean can't think of an easy inbuilt function. But one way could be to get the coordinates at due North and South from your center point and the height distance, and East and West using the width distance, for example. This should give you four coordinate points on the rectangle shap, then the combination of lats/lons should give you coordinate value for each corner of your shape. See here for Python: stackoverflow.com/questions/4530943/… – saph_top Feb 3 at 18:07
  • @saph_top - look above. Thanks again. – Dean Feb 3 at 23:13
0

Thank you @saph_top for your help.

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]

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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