Forgotten areas with BBoxSplitter

Apps / Services Python SDK
1

I’m trying to use BBoxSplitter in a download image script when final images will be bigger than 5000 pixels side. But it didn’t download all areas for my original big area. Sometimes it forgot half area or 90% of it…
I must be making a mistake but where? I really don’t understand.

Here is the part of my code who manage download for big areas:

blocks = gpd.read_file(blocks_name)

Xmin_m, Ymin_m, Xmax_m, Ymax_m = blocks.to_crs({'init': 'epsg:3857'}).total_bounds
print('Emprise :{0}'.format(blocks.to_crs({'init': 'epsg:4326'}).total_bounds))
print('Emprise :{0}'.format(blocks.to_crs({'init': 'epsg:3857'}).total_bounds))
largeur =  int(np.ceil((Xmax_m - Xmin_m)/10))
hauteur = int(np.ceil((Ymax_m - Ymin_m)/10))
print('\nArea pixel size: {0} x {1}'.format(largeur,hauteur))

if largeur > 5000 or hauteur > 5000:
    json = '{0}.json'.format(os.path.basename(blocks_path).replace('.shp', ''))
    if os.path.exists(json):
        os.remove(json)
    with open(json,'w') as f:
    f.write(blocks.to_json())
    if largeur > 5000:
        L = int(np.ceil(largeur/5000))
        print('%s cellules de large' % (L))
    if hauteur > 5000:
        H = int(np.ceil(hauteur/5000))
        print('%s cellules de haut' % (H))

geo_json = read_data(json)
area = shape(geo_json['features'][0]['geometry'])

bbox_splitter = BBoxSplitter([area], CRS.WGS84, (L, H), reduce_bbox_sizes=True)  # bounding box will be split into grid of LxH bounding boxes
bbox_list = bbox_splitter.get_bbox_list()

for emprise in bbox_list:
    dates_list = DLimages()

DLimages fonction juste download an image from a bbox.

Thanks in advance for your help

17
171114×897 27.6 KB

2

Hi Timothee,

I think the problem is that you are using coordinates in Popular Web Mercator projection (EPSG:3857) to calculate number of pixels in the area. This would only work if coordinates would be in UTM crs.

3

Ok. I do this to calculate image size in pixel knowing that a pixel is 10 meters apart. And so know if my future image is over 5000 pixels side. In a second time, it allows me to know in how many sub-bbox I have to split my main bbox.
So, how can I calculate this in UTM crs ??

4

Yes, in your code you use Xmin_m, Ymin_m, Xmax_m, Ymax_m in EPSG:3857 to calculate largeur and hauteur. But that is incorrect - the formula you are using would only work if Xmin_m, Ymin_m, Xmax_m, Ymax_m would be in UTM coordinate reference system.

Hence you have to project your data into UTM. As you might know there is not one but 128 different EPSG codes for UTM time zones. Hence you have to find the one that corresponds to your region of interest.

By the way, in sentinelhub.geo_utils there are some functions that might help you.

5

Thank you ! I’ll try it !