diff --git a/scripts/global_120x120km_UTM_tiling-grid.py b/scripts/global_120x120km_UTM_tiling-grid.py new file mode 100644 index 0000000..b692e37 --- /dev/null +++ b/scripts/global_120x120km_UTM_tiling-grid.py @@ -0,0 +1,120 @@ +""" +this script generates a global 120x120km tiling grid in UTM. + +""" + + +''' +STEPS +1. load the UTM zones gpkg +2. load the global land masses gpkg +3. loop over the 120 UTM zones + a. load the full 120x120 km grid depending on the hemisphere and change the epsg to the UTM zone + b. load a copy of the UTM zone polygon and warp to correct epsg + c. select the grid cells intersecting the UTM zone polygon and delete non-needed + d. load a copy of the land masses polygon and warp to correct epsg + e. select the grid cells intersecting the land masses and dlete non-needed ones + f. run over the remaining 120x120km tiles and add the naming as well as + add the openEO bbox dict using the bounds of the polygon (plus add the bounds in WGS84) + g. convert the geodataframe back to EPSG:4326 and save it to a new gpkg +4. combine all UTM zones into one gpkg and clean up the metadata +5. save the final global processing grid in the package ressources + +''' + +import geopandas as gpd +import pandas as pd +import math +from eo_processing.utils.mgrs import UTM_2_grid20id +import os + +# standard paths +path_grid_n = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\basic_120x120km_grid_nothern_no-crs_v2.gpkg') +path_grid_s = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\basic_120x120km_grid_southern_no-crs_v2.gpkg') +path_land = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\land_sea_mask_20kmbuffered_EPSG4326_v2.gpkg') +path_utm = os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\UTM_zones_high-res_EPSG4326.gpkg') + +# load UTM zones +gdf_utm = gpd.read_file(path_utm) + +# loop over all UTM zones +lZones = gdf_utm.name.unique().tolist() +lFiles = [] + +for UTMzone in lZones: + print(f'processing zone: {UTMzone}') + # get the UTM zone polygon and the EPSG + if UTMzone[-1] == 'N': + epsg = 32600 + int(UTMzone[:2]) + else: + epsg = 32700 + int(UTMzone[:2]) + + gdf_zone = gdf_utm[gdf_utm.name == UTMzone].copy() + + path_out = os.path.normpath( + r'C:\Users\buchhorm\Downloads\120x120km_grid\results\UTM_zone_{0}.gpkg'.format( + epsg)) + os.makedirs(os.path.dirname(path_out), exist_ok=True) + + if not os.path.exists(path_out): + # get the land masses + gdf_land = gpd.read_file(path_land) + gdf_land = gdf_land.clip(gdf_zone) + gdf_land = gdf_land.dissolve() + + # bring all to correct epsg + gdf_zone = gdf_zone.to_crs(epsg=epsg) + gdf_land = gdf_land.to_crs(epsg=epsg) + + # short cut if the gdf_land is empty + if gdf_land.empty: + continue + + # now we load the basic 120x120km grid and set the crs to the correct epsg + if UTMzone[-1] == 'N': + gdf_grid = gpd.read_file(path_grid_n).set_crs(epsg=epsg, allow_override=True) + else: + gdf_grid = gpd.read_file(path_grid_s).set_crs(epsg=epsg, allow_override=True) + + # select grid cells which are intersecting with utm zone + gdf_grid = gdf_grid[gdf_grid.intersects(gdf_zone.union_all(method='coverage'))] + + # intersecting + gdf_grid = gdf_grid[gdf_grid.intersects(gdf_land.union_all(method='coverage'))] + + # add the bbox dict for openEO + gdf_grid['bbox_dict'] = gdf_grid.apply( + lambda row: { + 'west': math.floor(row.geometry.bounds[0] / 100) * 100., + 'south': math.floor(row.geometry.bounds[1] / 100) * 100., + 'east': math.ceil(row.geometry.bounds[2] / 100) * 100., + 'north': math.ceil(row.geometry.bounds[3] / 100) * 100., + 'crs': epsg + }, axis=1 + ) + + # save to disk + + + gdf_grid[['left', 'top', 'right', 'bottom', 'row_index','col_index', 'bbox_dict', 'geometry']].to_file( + path_out) + else: + # load + gdf_grid = gpd.read_file(path_out) + gdf_zone = gdf_zone.to_crs(epsg=epsg) + gdf_grid = gdf_grid.clip(gdf_zone) + gdf_grid.to_file(path_out) + + lFiles.append(path_out) + + +# load all geopackages which were produced and convert to EPSG:4326 and merge them into one GeoDataFrame +result = gpd.GeoDataFrame(columns=['left', 'top', 'right', 'bottom', 'row_index','col_index', 'bbox_dict', 'geometry'], geometry='geometry', crs='EPSG:4326') + +for file in lFiles: + gdf_tmp = gpd.read_file(file) + gdf_tmp = gdf_tmp.to_crs(epsg=4326) + result = pd.concat([result, gdf_tmp], ignore_index=True) + +# save final result to disk +result.to_file(os.path.normpath(r'C:\Users\buchhorm\Downloads\120x120km_grid\ECDC_global_120x120km_grid.gpkg')) diff --git a/scripts/laea_tiling_grids_2_high_res_Sentinel2_tiles.py b/scripts/laea_tiling_grids_2_high_res_Sentinel2_tiles.py index eb6500a..4444e26 100644 --- a/scripts/laea_tiling_grids_2_high_res_Sentinel2_tiles.py +++ b/scripts/laea_tiling_grids_2_high_res_Sentinel2_tiles.py @@ -28,17 +28,17 @@ def get_epsg(tile_id): return target_EPSG # filter to tiles in UTM zones needed for pan EU -df2 = pd.read_csv(r'C:\Users\buchhorm\Downloads\new_grids\utm_codes_europe.csv') -df2['utm'] = df2['ZONE'].astype(str) + df2['ROW_'] -lNeeded = df2['utm'].tolist() +#df2 = pd.read_csv(r'C:\Users\buchhorm\Downloads\new_grids\utm_codes_europe.csv') +#df2['utm'] = df2['ZONE'].astype(str) + df2['ROW_'] +#lNeeded = df2['utm'].tolist() ## Start read out everything we need with beautiful soup from Sentinel-2 KML -xml_file = r'C:\Users\buchhorm\Downloads\new_grids\Sentinel2_tiles.kml' -soup = bs.BeautifulSoup(open(xml_file), 'xml') +xml_file = r'C:\Users\buchhorm\Downloads\120x120km_grid\S2_tiles.kml' +soup = bs.BeautifulSoup(open(xml_file), 'html.parser') #find all VRTRasterBands in the soup -tiles = soup.findAll('description') +tiles = soup.find_all('description') data = [] @@ -57,8 +57,8 @@ def get_epsg(tile_id): cols = [ele.text.strip() for ele in cols] if cols[0] == 'TILE_ID': data1 = cols[1] if cols[0] == 'UTM_WKT': data2 = cols[1] - if data1[:3] in lNeeded: - data.append((data1, data2, get_epsg(data1))) + + data.append((data1, data2, get_epsg(data1))) df = pd.DataFrame(data, columns=['tile_id', 'geometry', 'epsg']) data = None @@ -80,10 +80,10 @@ def get_epsg(tile_id): # add extra points in 250m intervall gdf['geometry'] = gdf['geometry'].apply(lambda x: x.segmentize(250)) - gdf2 = gdf.to_crs(epsg=3035) + gdf2 = gdf.to_crs(epsg=4326) list_results.append(gdf2) result = gpd.GeoDataFrame(pd.concat(list_results, ignore_index=True), crs=list_results[0].crs) # write out -result.to_file(r'C:\Users\buchhorm\Downloads\new_grids\Sentinel2_tiling_grid_EU_high_res_EPSG3035.gpkg', driver='GPKG') +result.to_file(r'C:\Users\buchhorm\Downloads\120x120km_grid\Sentinel2_tiling_grid_high_res_EPSG4326.gpkg', driver='GPKG') diff --git a/src/eo_processing/resources/ECDC_global_120x120km_grid_v1.gpkg b/src/eo_processing/resources/ECDC_global_120x120km_grid_v1.gpkg new file mode 100644 index 0000000..45f7762 Binary files /dev/null and b/src/eo_processing/resources/ECDC_global_120x120km_grid_v1.gpkg differ