Max distance between polygons in a polygon set

I have a polygon set in GIS with a lot of polygons. My task is to find two polygons with max distance between them [two polygons that are the farthest/farthermost, my english is funny :( ]

This is my code so far, but I stopped because output gives me the the max value that is not correct. Help please. This is my code:

# -*- coding: utf-8 -*-
from shapely.wkb import loads
from osgeo import ogr


datoteka=ogr.Open("sth.shp")
sloj=datoteka.GetLayerByName("katastar")


udaljenosti=[]
for i in sloj:
    cestica1=loads(i.GetGeometryRef().ExportToWkb())
    for j in sloj:
        cestica2=loads(j.GetGeometryRef().ExportToWkb())
        udaljenosti.append(cestica1.distance(cestica2))


max_udaljenost=max(udaljenosti)
print max_udaljenost


datoteka.Destroy()

i know i have a lot of unnecessary stuff in the code but I'll fix that when i fix this maximum distance.

Solution

There is most likely an issue with how you iterate over the layer sloj. The thing is that if you write a code like this:

#!/usr/bin/env python
from shapely.wkb import loads
from osgeo import ogr

datoteka = ogr.Open("custom.shp")
sloj = datoteka.GetLayer()

N = sloj.GetFeatureCount()
print(N)
for p in sloj:
    cnt = 0
    for q in sloj:
        cnt += 1
    print(cnt)

then this won't print N times the number N, but instead only number N followed by N-1. The reason is that sloj works as an iterator, thus when the outer loop is entered, it consumes the first element and the inner loop consumes all remaining N-1 elements. Therefore the outer loop terminates on the next iteration since there is nothing left to iterate over.

In order to reflect this in your code, you could for example do (if the number of objects is not too large):

# -*- coding: utf-8 -*-
from shapely.wkb import loads
from osgeo import ogr


datoteka=ogr.Open("C:\Users\Sanja\Desktop\Desktop\Geof\Geof 
IV\BPP\podaci\katastar.shp")
sloj=datoteka.GetLayerByName("katastar")

L = [p.GetGeometryRef().ExportToWkb() for p in sloj]
N, max_dist = len(L), -1

for i in range(N):
    cestica1 = L[i]
    for j in range(i+1, N):
        cestica2 = L[j]
        max_dist = max(max_dist, cestica1.distance(cestica2))

print(max_dist)

datoteka.Destroy()

However, these distances are calculated in the lat/lon coordinates, i.e., they do not represent true distances on the globe. Nevertheless, if the objects in your layer don't span large areas of the globe, then the pair with the largest distance in the lat/lon coordinate system will most likely correspond to the pair with the largest "true" distance.