2

I am using the EE API in a Python Jupyter Notebook.

I am building a map based on the LSIB_SIMPLE country boundaries. My objective is to display my information in a map and adapt the zoom to fit the country size (Singapore is way smaller than Congo for example)

I would like to compute both the longest north-south and east-west lines inside the country to adapt my zoom accordingly on this object:

country_code = CG
country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017').filter(ee.Filter.eq('country_co', country_code))

geom = country.geometry()
  • Is there a way to find the 4 needed points in the geometry object (northest, southest, eastest, westest) so that I can compute distances between them ?
6
  • what does the notebook show after you type geom. and then hit the tab key?
    – Paul H
    Jul 7, 2020 at 22:17
  • nothing actually Jul 7, 2020 at 22:21
  • how about dir(geom)
    – Paul H
    Jul 7, 2020 at 22:29
  • too many thing to go in a comment... I checked .coordinates that is an invocation but with the .getInfo()["coordinates"] I retreived the array of coordinates. Jul 7, 2020 at 22:36
  • you could throw those into a shapely geometry and get the area and extent
    – Paul H
    Jul 8, 2020 at 1:08

3 Answers 3

1

thanks to @Paul H for the idea.

I ended up finding a solution that compute the 4 cardinal points of my AOI and adapt the zoom to the longest diagonal :


from haversine import haversine

def update_zoom(asset_id):
    """search for the dimension of the AOI and adapt the map zoom acordingly
    
    Args: 
        asset_id (str): the assetID
    
    Returns: 
        zoom (int): the zoom value riquired
    """
    
    #retreive the asset 
    geom = ee.FeatureCollection(asset_id).geometry()
    coordinates = geom.getInfo()["coordinates"]
    #transform into a single list of all the coordinates
    shape = []
    get_coords(coordinates, shape)
    #in the coordinates search for the 4 cardinal points of the aoi
    #gee format coords [lng, lat]
    count = 0 
    for coords in shape: #perimeter of each shape
        count += 1
        if count == 1:
            north = coords
            east = coords
            south = coords
            west = coords
            continue
            
        if coords[0] < west[0]:
            west = coords
        if coords[0] > east[0]:
            east = coords
        if coords[1] < south[1]:
            south = coords
        if coords[1] > north[1]:
            north = coords

    maxsize = max(haversine(east, west), haversine(north, south))
    
    lg = 40075 #number of displayed km at zoom 1
    zoom = 1
    while lg > maxsize:
        zoom += 1
        lg /= 2
        
    return zoom-1

def get_coords(coordinates, array_coord=[]):
    """get all the coordinates and set them in a single table of tuple without knowing in advance the depth of the tab"""
    if isinstance(coordinates[0], float):
        array_coord.append((coordinates[0], coordinates[1]))
        return array_coord
    
    for item in coordinates:
        get_coords(item,array_coord)

feel free to suggest modifications there is always room for improvement

2
1

You can find min/max coordinates easily with server-side Earth Engine. Call .bounds() on the geometry or feature, get the .coordinates(), then reduce the lats and lons by min and max.

geometry = ee.Geometry.Polygon(
        [[[-121.04086404798966, 36.97562343255799],
          [-121.38144021986466, 36.93611720712966],
          [-121.06283670423966, 36.544342687193584],
          [-120.47506814955216, 36.38088239731507],
          [-120.10153299330216, 36.50902912598594],
          [-120.46957498548966, 36.486949964934816],
          [-120.52999979017716, 36.96245696791018],
          [-119.93124490736466, 37.14658010572364],
          [-120.62887674330216, 37.51348410687124],
          [-121.16171365736466, 37.16409233412239]]]);

bounds = ee.Array(ee.List(geometry.bounds().coordinates()).get(0))

min_coords = bounds.reduce(ee.Reducer.min(), [0]).project([1]).toList()
max_coords = bounds.reduce(ee.Reducer.max(), [0]).project([1]).toList()

x_min = min_coords.get(0)
y_min = min_coords.get(1)
x_max = max_coords.get(0)
y_max = max_coords.get(1)

from pprint import pprint
pprint(bounds.getInfo())
print()
print('x_min:', x_min.getInfo())
print('y_min:', y_min.getInfo())
print('x_max:', x_max.getInfo())
print('y_max:', y_max.getInfo())
0

With Python the script can be as succinct as below:

def MinMaxXY(geometry):
    x_min = ee.List(ee.List(geometry.coordinates().get(0)).map(lambda x: ee.List(x).get(0))).reduce(reducer=ee.Reducer.min()).getInfo()
    x_max = ee.List(ee.List(geometry.coordinates().get(0)).map(lambda x: ee.List(x).get(0))).reduce(reducer=ee.Reducer.max()).getInfo()
    y_min = ee.List(ee.List(geometry.coordinates().get(0)).map(lambda x: ee.List(x).get(1))).reduce(reducer=ee.Reducer.min()).getInfo()
    y_max = ee.List(ee.List(geometry.coordinates().get(0)).map(lambda x: ee.List(x).get(1))).reduce(reducer=ee.Reducer.max()).getInfo()
    return x_min, x_max, y_min, y_max

And a simple test:

geometry = ee.Geometry.Polygon(
        [[[-121.04086404798966, 36.97562343255799],
          [-121.38144021986466, 36.93611720712966],
          [-121.06283670423966, 36.544342687193584],
          [-120.47506814955216, 36.38088239731507],
          [-120.10153299330216, 36.50902912598594],
          [-120.46957498548966, 36.486949964934816],
          [-120.52999979017716, 36.96245696791018],
          [-119.93124490736466, 37.14658010572364],
          [-120.62887674330216, 37.51348410687124],
          [-121.16171365736466, 37.16409233412239]]])

MinMaxXY(geometry)

(-121.38144021986466,
 -119.93124490736466,
 36.38088239731507,
 37.51348410687124)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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