Creating 1 hectare Polygons given latitude longitude in Python

Question

I'm trying to work with openweathermap's Agromonitoring API.

To use their API you have to create polygons representing your area, then you can request weather data, soil moisture etc for that particular area/polygon. There are 2 ways of creating polygons:

1. Use their graphical tool to draw the polygon on a map manually.
2. Use their polygon API to create the polygon given a list of coordinates.

I'm interested in the second method, simple get requests with the following body accomplish polygon creation for me.

body = {
"name":"TestPolygon",
"geo_json":{
   "type":"Feature",
   "properties":{

   },
   "geometry":{
      "type":"Polygon",
      "coordinates":[
         [[0, 0], [0, 0.0009], [0.0009, 0.0009], [0.0009, 0], [0, 0]]
      ]
   }
  }
}

Their API requires a minimum polygon area of 1 hectare, so I tried this request with the above dummy coordinates and the polygon created had an area of 1.003 hectares. I want it to be as close to one as possible.

How do I create an appropriate coordinates list given a single point (lat, lng), which would result in a polygon of 1 hectare?

I only want to input latitude and longitude, and then get a polygon as close to 1 hectare as possible. Perhaps the single point can be the center of our polygon (and we expand around it), or it can be a corner point of a rectangle shape. The shape doesn't matter exactly, it just needs to be nearly 1 hectare.

Since, my coordinates list brought me a 1 hectare polygon with 0.0009, I thought I could fix 0.0009 to be my factor for all coordinates, like so:

# I tried creating a rectangular polygon.
lat = 25.46666
lng = 68.71666
factor = 0.0009
lat_lng_list = [[lng, lat], [lng, lat + factor], [lng + factor, lat + factor], [lng + factor, lat], [lng, lat]]

but to my surprise, this factor does not remain constant for all latitude,longitude points, And it needs to be changed (through trial and error) until I reach a 1 hectare polygon for that specific point. It's also worth noting that the only way I find out if my polygon is 1 hectare or close to it, is through the response I receive to my get request.

GIS is not my field. I'm interested in hearing how you would solve this problem or how I can use code to automate polygon creation (so the user won't have to manually draw it on their map).

Just to reiterate the problem statement: Given a single point (lat, lng), create a 1 hectare polygon of any shape with that point as the center or as the corner (whichever is easier I have no clue).

Answer 1

Here's an idea, creating a function which automatically creates a polygon with the area of 1 hectare for you according to the following steps:

1. Input the central coordinate
2. Reproject to an equal-area projection
3. Create a Square with the area of exactly 1 hectare
4. Project back to WGS 84

Here's an example implementation outputting a GeoJSON containing the desired polygon:

import geopandas as gpd
def create_hectar_sized_square(lat, lon):
    gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(lon, lat)])
    gdf.crs = {'init': 'epsg:4326'}
    gdf = gdf.to_crs({'init': 'epsg:5070'})
    gdf['geometry'] = gdf['geometry'].apply(lambda x: x.buffer(50).envelope)
    print(gdf.area)
    gdf = gdf.to_crs({'init': 'epsg:4326'})
    gdf.plot()
    return gdf.to_json()

An here is an example usage:

create_hectar_sized_square(53.43079, -2.96071)

which will output:

'{"type": "FeatureCollection",
 "features": [{"id": "0",
               "type": "Feature",
               "properties": {},
               "geometry": {"type": "Polygon",
                            "coordinates": [
                                [[-2.9611591576420615, 53.43052239401445], 
                                 [-2.9602608423579424, 53.43052239401445],
                                 [-2.9602608423579424, 53.4310576043007],
                                 [-2.9611591576420615, 53.4310576043007],
                                 [-2.9611591576420615, 53.43052239401445]]
                             ]}
                }]
   }'

Edit: As noted by Keith Ma, it's important to chose an appropriate equal-area projection for your area of interest. Answer updated accordingly.

Answer 2

Note that bexi's solution picks an unfortunate coordinate system to perform the buffer operation in. Web mercator is not an equal area projection, and so the area of the reprojected geometry is not one hectare. The difference is not small either. For a central point in Ohio, USA, I get an area of ~0.58 ha.

If you replace epsg:3857 with a suitable equal area projection (say epsg:5070 for the continental US), the method works like a charm.