3

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.

Now the actual question, is 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, so I apologize if this question is silly, but 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).

  • The Earth is an oblate spheroid. All calculations on that spheroid will, of necessity, be somewhat rough. You'd need to use an infinite number of vertices to obtain a perfect circle; using a reasonable number of vertices will introduce error. I'm working a project right now that takes 3-20 iterations to achieve "five nines" accuracy on 25 sq. km "rectangles" (geodetic trapezoids), and I have a large geodetic toolbox at my disposal. There's no reason you can't script your current manual process. If you want more that that, you may need to immerse yourself in geodesy. – Vince Aug 6 at 12:12
  • @Vince I did consider scripting my current manual process, but wanted to know if there's possibly a better/simpler way to achieve this. I'm not looking for perfect accuracy or perfect shapes so no worries, I imagine there must be a lot of approximations in these types of calculations anyway. I'll google geodesy just to find out what it means haha, thank you for your input! – AbdurRehman Khan Aug 6 at 12:19
0

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 a projection using metres as units (e.g. WebMercator)
  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:3857'})
    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]]
                             ]}
                }]
   }'
  • Wow, this is quite brilliant! Thank you, I learned a lot from your answer :) Definitely motivated to explore geopandas now. – AbdurRehman Khan Aug 19 at 11:43

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.