diff --git a/data/KaartbladenNGI/Kbl.lyr b/data/KaartbladenNGI/Kbl.lyr new file mode 100644 index 0000000..470543d Binary files /dev/null and b/data/KaartbladenNGI/Kbl.lyr differ diff --git a/data/KaartbladenNGI/Kbl.prj b/data/KaartbladenNGI/Kbl.prj new file mode 100644 index 0000000..2c385eb --- /dev/null +++ b/data/KaartbladenNGI/Kbl.prj @@ -0,0 +1 @@ +PROJCS["Belge_Lambert_1972",GEOGCS["GCS_Belge_1972",DATUM["D_Belge_1972",SPHEROID["International_1924",6378388.0,297.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",150000.01256],PARAMETER["False_Northing",5400088.4378],PARAMETER["Central_Meridian",4.367486666666666],PARAMETER["Standard_Parallel_1",49.8333339],PARAMETER["Standard_Parallel_2",51.16666723333333],PARAMETER["Latitude_Of_Origin",90.0],UNIT["Meter",1.0]] \ No newline at end of file diff --git a/source/Orthophoto_processing.Rmd b/source/Orthophoto_processing.Rmd index 954483e..36f69dc 100644 --- a/source/Orthophoto_processing.Rmd +++ b/source/Orthophoto_processing.Rmd @@ -3,60 +3,109 @@ title: "Orthophoto analysis" output: html_notebook --- -```{r get orthophotos through WMTS} -library(xml2) +```{r match drone plots with the correct orthophoto tiles} +library(terra) +library(sf) -# Your WMTS GetCapabilities URL -wmts_url <- "https://geo.api.vlaanderen.be/omz/wmts?request=getcapabilities&service=wmts&version=1.0.0" +# Load the GeoJSON file of the drone plot outlines +drone_outlines <- st_read("../data/plots.geojson") +print(drone_outlines) -# Fetch the WMTS capabilities XML -wmts_capabilities <- read_xml(wmts_url) +# Load the tiling shapefile (kaartbladen NGI) +tile_shapefile <- st_read("../data/KaartbladenNGI/Kbl.shp") -# Parse information from the XML -subnodes <- xml_children(wmts_capabilities) +# Add the "VL_" prefix to the CODE column +tile_shapefile$TileName <- paste0("OMZCIR21VL_", sprintf("%02d", as.integer(tile_shapefile$CODE))) -# get into the Contents part of the XML -contents_node <- subnodes[[4]] +# Ensure both datasets have the same CRS +drone_outlines <- st_transform(drone_outlines, st_crs(tile_shapefile)) -# Get all nodes under -layers <- xml_children(contents_node) +# Spatial join to find intersections +matched_tiles <- st_join(drone_outlines, tile_shapefile, left = FALSE) -# Loop through each and extract information -layers_info <- lapply(layers, function(layer_node) { - title <- xml_text(xml_find_first(layer_node, ".//ows:Title")) - identifier <- xml_text(xml_find_first(layer_node, ".//ows:Identifier")) +# Inspect the result +print(matched_tiles) + + +# General folder path for orthos (GIS schijf INBO) +ortho_folder <- "Z:/Vlaanderen/Referentie/Orthofoto/Orthofoto_zomer2021_CIR" +## List ortho files (tiles) +ortho_files <- list.files(ortho_folder, pattern = "\\.jp2$", full.names = TRUE) + +# Extract tile names from raster filenames +ortho_data <- data.frame( + TileName = gsub("\\.jp2", "\\1", basename(ortho_files)), + RasterFile = ortho_files +) + +# Merge raster file paths with the matched tiles +matched_tiles <- merge(matched_tiles, ortho_data, by = "TileName", all.x = TRUE) + +print(matched_tiles) +``` +```{r extract orthophoto parts overlapping the drone plots} +# Loop through each unique drone polygon +for (drone_id in unique(matched_tiles$naam)) { + # Subset the current drone polygon + drone_geom <- matched_tiles[matched_tiles$naam == drone_id, ] + print(drone_id) - # Bounding box - bbox_lower <- xml_text(xml_find_first(layer_node, ".//ows:LowerCorner")) - bbox_upper <- xml_text(xml_find_first(layer_node, ".//ows:UpperCorner")) + # Add a 10-meter buffer around the polygon + drone_geom <- st_buffer(drone_geom, dist = 10) - # Format - format <- xml_text(xml_find_first(layer_node, ".//*[local-name()='Format']")) - # Tile matrix sets - tile_matrix_sets <- paste(xml_text(xml_find_all(layer_node, ".//*[local-name()='TileMatrixSet']")), collapse = ", ") + # Find all raster tiles corresponding to this drone polygon + tile_files <- unique(drone_geom$RasterFile) + (print(tile_files)) + # Initialize a list to store clipped rasters + clipped_parts <- list() - # Return as a list - data.frame( - Title = title, - Identifier = identifier, - BoundingBoxLower = bbox_lower, - BoundingBoxUpper = bbox_upper, - Format = format, - TileMatrixSets = tile_matrix_sets, - stringsAsFactors = FALSE - ) -}) - -# Combine all layer data frames into one -layers_df <- do.call(rbind, layers_info) - -# Print the resulting data frame -print(layers_df) - -``` -```{r get forest polygons from WMTS} - - + # Loop through the corresponding raster tiles + for (tile_file in tile_files) { + # Check if the raster file exists + if (!file.exists(tile_file)) { + cat("Raster file not found:", tile_file, "\n") + next + } + + # Load the raster + raster <- rast(tile_file) + + # Clip and mask the raster to the drone polygon + print('starting clip') + # Reproject the vector to match the raster + crs_raster <- crs(raster) + drone_geom <- st_transform(drone_geom, crs_raster) + clipped_raster <- crop(raster, vect(drone_geom)) + masked_raster <- mask(clipped_raster, vect(drone_geom)) + + # Add to the list of clipped parts + clipped_parts <- append(clipped_parts, list(masked_raster)) + } + print('clipping done') + print('starting merge') + # Merge all clipped parts if there are multiple tiles + if (length(clipped_parts) > 1) { + merged_raster <- do.call(merge, clipped_parts) + } else if (length(clipped_parts) == 1) { + merged_raster <- clipped_parts[[1]] + } else { + cat("No valid raster data for Drone ID:", drone_id, "\n") + next + } + print('merge done') + # Define the output file name + output_file <- paste0("../output/ortho_CIR/ortho_CIR_", drone_id, "_10mbuf.tif") + # Ensure the output directory exists + if (!dir.exists("../output/ortho_CIR")) { + dir.create("../output/ortho_CIR", recursive = TRUE) + } + print('starting save raster') + # Save the merged raster + writeRaster(merged_raster, output_file, overwrite = TRUE) + print('save raster done') + # Print status + cat("Saved merged raster for Drone ID:", drone_id, "to", output_file, "\n") +} ```