@@ -267,53 +267,63 @@ def convert_geotiff_to_geotiff(
267
267
nodata_value ,
268
268
)
269
269
270
- def generate_data (self ) -> np .ndarray :
271
- """Generate data from the NRW provider.
270
+ def transform_bbox (self , bbox : tuple [float , float , float , float ], to_crs : str ) -> tuple [float , float , float , float ]:
271
+ """Transform the bounding box to a different coordinate reference system (CRS).
272
+
273
+ Arguments:
274
+ bbox (tuple): Bounding box to transform (north, south, east, west).
275
+ to_crs (str): Target CRS (e.g., EPSG:4326 for CRS:84).
272
276
273
277
Returns:
274
- np.ndarray: Numpy array of the data .
278
+ tuple: Transformed bounding box (north, south, east, west) .
275
279
"""
276
- wcs = WebCoverageService ('https://www.wcs.nrw.de/geobasis/wcs_nw_dgm' , auth = Authentication (verify = False ), timeout = 600 )
277
- north , south , east , west = self .get_bbox ()
278
-
279
- # Because the target CRS is rotated compared to 4326, we need to transform all corners and get the bounding box that covers the rotated area
280
- transformer = Transformer .from_crs ("epsg:4326" , "epsg:25832" )
280
+ transformer = Transformer .from_crs ("epsg:4326" , to_crs )
281
+ north , south , east , west = bbox
281
282
bottom_left_x , bottom_left_y = transformer .transform (xx = south , yy = west )
282
283
top_left_x , top_left_y = transformer .transform (xx = north , yy = west )
283
284
top_right_x , top_right_y = transformer .transform (xx = north , yy = east )
284
285
bottom_right_x , bottom_right_y = transformer .transform (xx = south , yy = east )
285
286
286
- # get the bounding box that covers the rotated area
287
- # coordinate system is reversed, so this looks weird but it is correct
288
287
west = min (bottom_left_y , bottom_right_y )
289
288
east = max (top_left_y , top_right_y )
290
289
south = min (bottom_left_x , top_left_x )
291
290
north = max (bottom_right_x , top_right_x )
292
291
293
- print ( west , east , south , north )
292
+ return north , south , east , west
294
293
295
- # Generate coordinates for the tile corners
296
- x_coords = np .arange (west , east , 1000 )
297
- y_coords = np .arange (south , north , 1000 )
294
+ def tile_bbox (self , bbox : tuple [float , float , float , float ], tile_size : int ) -> list [tuple [float , float , float , float ]]:
295
+ """Tile the bounding box into smaller bounding boxes of a specified size.
298
296
299
- # Append the edge coordinates because arange does not include the end value
297
+ Arguments:
298
+ bbox (tuple): Bounding box to tile (north, south, east, west).
299
+ tile_size (int): Size of the tiles in meters.
300
+
301
+ Returns:
302
+ list: List of smaller bounding boxes (north, south, east, west).
303
+ """
304
+ north , south , east , west = bbox
305
+ x_coords = np .arange (west , east , tile_size )
306
+ y_coords = np .arange (south , north , tile_size )
300
307
x_coords = np .append (x_coords , east ).astype (x_coords .dtype )
301
308
y_coords = np .append (y_coords , north ).astype (y_coords .dtype )
302
309
303
- print (x_coords )
304
- print (y_coords )
305
-
306
- # Use meshgrid to create grid combinations
307
310
x_min , y_min = np .meshgrid (x_coords [:- 1 ], y_coords [:- 1 ], indexing = "ij" )
308
311
x_max , y_max = np .meshgrid (x_coords [1 :], y_coords [1 :], indexing = "ij" )
309
312
310
- # Combine into a single array of bounding boxes
311
313
tiles = np .stack ([x_min .ravel (), y_min .ravel (), x_max .ravel (), y_max .ravel ()], axis = 1 )
312
314
313
- print ( tiles )
315
+ return tiles
314
316
315
- all_tif_files = []
317
+ def download_tiles (self , tiles : list [tuple [float , float , float , float ]]) -> list [str ]:
318
+ """Download tiles from the NRW provider.
319
+
320
+ Arguments:
321
+ tiles (list): List of tiles to download.
316
322
323
+ Returns:
324
+ list: List of paths to the downloaded GeoTIFF files.
325
+ """
326
+ all_tif_files = []
317
327
for tile in tiles :
318
328
file_name = '_' .join (map (str , tile )) + '.tif'
319
329
file_path = os .path .join (self .shared_tiff_path , file_name )
@@ -328,6 +338,21 @@ def generate_data(self) -> np.ndarray:
328
338
f .write (output .read ())
329
339
330
340
all_tif_files .append (file_path )
341
+ return all_tif_files
342
+
343
+ def generate_data (self ) -> np .ndarray :
344
+ """Generate data from the NRW provider.
345
+
346
+ Returns:
347
+ np.ndarray: Numpy array of the data.
348
+ """
349
+ north , south , east , west = self .get_bbox ()
350
+
351
+ north , south , east , west = self .transform_bbox ((north , south , east , west ), "epsg:25832" )
352
+
353
+ tiles = self .tile_bbox ((north , south , east , west ), 1000 )
354
+
355
+ all_tif_files = self .download_tiles (tiles )
331
356
332
357
self .merge_geotiff (all_tif_files , os .path .join (self .output_path , "merged.tif" ))
333
358
0 commit comments