Skip to content
Snippets Groups Projects
Commit b875c9e3 authored by Zvolenský Juraj's avatar Zvolenský Juraj
Browse files

ran formatter

parent 572dae8d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
swe_maps_stac/compare_data.ipynb
+ 63
− 39
View file @ b875c9e3
%% Cell type:markdown id: tags:
## Comparing the data before and after
For quality control, the original NetCDF data is compared to the newly created COGs. \
This is to ensure that the spatial referencing and data values are consistent between the two formats. \
%% Cell type:code id: tags:
``` python
import os
import xarray as xr
import numpy as np
import rasterio
import matplotlib.pyplot as plt
from datetime import datetime
```
%% Cell type:code id: tags:
``` python
cog_folders = [
'./output_cog/venosta_1718_cog',
'./output_cog/venosta_1819_cog',
'./output_cog/venosta_1920_cog',
'./output_cog/venosta_2021_cog',
'./output_cog/venosta_2122_cog',
"./output_cog/venosta_1718_cog",
"./output_cog/venosta_1819_cog",
"./output_cog/venosta_1920_cog",
"./output_cog/venosta_2021_cog",
"./output_cog/venosta_2122_cog",
]
sample_cogs = [
'./output_cog/venosta_1718_cog/swe_reconstruction_c_50m_s_20171002_it_epsg.3035_v20240212',
'./output_cog/venosta_1819_cog/swe_reconstruction_c_50m_s_20180121_it_epsg.3035_v20240212',
'./output_cog/venosta_1920_cog/swe_reconstruction_c_50m_s_20191201_it_epsg.3035_v20240212',
'./output_cog/venosta_2021_cog/swe_reconstruction_c_50m_s_20200515_it_epsg.3035_v20240212',
'./output_cog/venosta_2122_cog/swe_reconstruction_c_50m_s_20210605_it_epsg.3035_v20240212',
"./output_cog/venosta_1718_cog/swe_reconstruction_c_50m_s_20171002_it_epsg.3035_v20240212.tif",
"./output_cog/venosta_1819_cog/swe_reconstruction_c_50m_s_20181005_it_epsg.3035_v20240212.tif",
"./output_cog/venosta_1920_cog/swe_reconstruction_c_50m_s_20191201_it_epsg.3035_v20240212.tif",
"./output_cog/venosta_2021_cog/swe_reconstruction_c_50m_s_20210215_it_epsg.3035_v20240212.tif",
"./output_cog/venosta_2122_cog/swe_reconstruction_c_50m_s_20220412_it_epsg.3035_v20240212.tif",
]
```
%% Cell type:code id: tags:
``` python
def count_files_in_dirs(dir_paths: list) -> None:
for dir_path in dir_paths:
if os.path.isdir(dir_path):
num_files = len([f for f in os.listdir(dir_path) if os.path.isfile(os.path.join(dir_path, f))])
num_files = len(
[
f
for f in os.listdir(dir_path)
if os.path.isfile(os.path.join(dir_path, f))
]
)
print(f"There are {num_files} files in {dir_path}")
else:
print(f"{dir_path} is not a valid directory.")
count_files_in_dirs(cog_folders)
```
%% Output
There are 365 files in ./output_cog/venosta_1718_cog
There are 365 files in ./output_cog/venosta_1819_cog
There are 366 files in ./output_cog/venosta_1920_cog
There are 365 files in ./output_cog/venosta_2021_cog
There are 365 files in ./output_cog/venosta_2122_cog
%% Cell type:code id: tags:
``` python
def print_cog_info(cog_paths: list) -> None:
for cog_path in cog_paths:
if os.path.isfile(cog_path):
with rasterio.open(cog_path) as src:
print(f"File: {cog_path}")
print(f"CRS: {src.crs}")
print(f"Size: {src.width} x {src.height}")
print(f"Bounds: {src.bounds}")
# Extract date from filename, assuming it's in the format "YYYYMMDD"
filename = os.path.basename(cog_path)
date_str = ''.join(filter(str.isdigit, filename))
if len(date_str) >= 8:
date = datetime.strptime(date_str[:8], "%Y%m%d")
print(f"Date: {date}")
else:
print(f"{cog_path} is not a valid file.")
print_cog_info(sample_cogs)
```
%% Output
File: ./output_cog/venosta_1718_cog/swe_reconstruction_c_50m_s_20171002_it_epsg.3035_v20240212.tif
CRS: EPSG:3035
Size: 1233 x 974
Bounds: BoundingBox(left=4344231.07781725, bottom=2592093.7705513844, right=4405838.178041091, top=2640759.882163729)
File: ./output_cog/venosta_1819_cog/swe_reconstruction_c_50m_s_20181005_it_epsg.3035_v20240212.tif
CRS: EPSG:3035
Size: 1233 x 974
Bounds: BoundingBox(left=4344231.07781725, bottom=2592093.7705513844, right=4405838.178041091, top=2640759.882163729)
File: ./output_cog/venosta_1920_cog/swe_reconstruction_c_50m_s_20191201_it_epsg.3035_v20240212.tif
CRS: EPSG:3035
Size: 1233 x 974
Bounds: BoundingBox(left=4344231.07781725, bottom=2592093.7705513844, right=4405838.178041091, top=2640759.882163729)
File: ./output_cog/venosta_2021_cog/swe_reconstruction_c_50m_s_20210215_it_epsg.3035_v20240212.tif
CRS: EPSG:3035
Size: 1233 x 974
Bounds: BoundingBox(left=4344231.07781725, bottom=2592093.7705513844, right=4405838.178041091, top=2640759.882163729)
File: ./output_cog/venosta_2122_cog/swe_reconstruction_c_50m_s_20220412_it_epsg.3035_v20240212.tif
CRS: EPSG:3035
Size: 1233 x 974
Bounds: BoundingBox(left=4344231.07781725, bottom=2592093.7705513844, right=4405838.178041091, top=2640759.882163729)
%% Cell type:code id: tags:
``` python
```
......
swe_maps_stac/conversion.py
+ 2
− 4
View file @ b875c9e3
......@@ -38,7 +38,7 @@ def netcdf_to_cog(netcdf_path: str, cog_dir: str) -> None:
dataset.sizes["y"],
)
for i in tqdm(range(dataset.sizes["time"])):
for i in tqdm(range(dataset.sizes["time"])):
time_slice = dataset.isel(time=i)
timestamp = time_slice.time.values.astype("datetime64[D]").astype(str)
......@@ -95,7 +95,5 @@ def netcdf_to_cog(netcdf_path: str, cog_dir: str) -> None:
def process_multiple_datasets(netcdf_paths, cog_dirs):
for netcdf_path, cog_dir in zip(
netcdf_paths, cog_dirs
):
for netcdf_path, cog_dir in zip(netcdf_paths, cog_dirs):
netcdf_to_cog(netcdf_path, cog_dir)
swe_maps_stac/convert.ipynb
+ 29
− 43
View file @ b875c9e3
%% Cell type:markdown id: tags:
## Converting NetCDF into COGs
Copying the data from the network drive using a wildcard to get SWE maps. \
Then using a function defined in utils.py to convert the NetCDF files into COGs.
%% Cell type:code id: tags:
``` python
from utils import copy_netcdf_data
from conversion import netcdf_to_cog, process_multiple_datasets
from transformation import transform_cogs_in_dirs
```
%% Cell type:code id: tags:
``` python
# Find and copy the SWE data to the data folder
netcdf_path = '/mnt/CEPH_PROJECTS/PROSNOW/SENTINEL-2/Daily_SCA/Venosta/*_swe.nc'
output = './data'
netcdf_path = "/mnt/CEPH_PROJECTS/PROSNOW/SENTINEL-2/Daily_SCA/Venosta/*_swe.nc"
output = "./data"
copy_netcdf_data(netcdf_path, output)
```
%% Cell type:code id: tags:
``` python
path = '/mnt/CEPH_PROJECTS/PROSNOW/SENTINEL-2/Daily_SCA/Venosta/venosta_hy2021_swe.nc'
output = './data'
copy_netcdf_data(path, output)
```
%% Cell type:code id: tags:
``` python
cogs= netcdf_to_cog(netcdf_path = './data/venosta_hy2021_swe.nc',
cog_dir = './output_cog/venosta_2021_cog/',
)
```
%% Cell type:code id: tags:
``` python
# Convert single NetCDF into a folder of COGs
cogs = netcdf_to_cog(
netcdf_path="./data/venosta_hy2021_swe.nc",
cog_dir="./output_cog/venosta_2021_cog/",
)
```
%% Cell type:code id: tags:
``` python
# Convert a list of NetCDF files into their respective output COG folders
net_cdf_paths = [
'./data/venosta_hy1718_swe.nc',
'./data/venosta_hy1819_swe.nc',
'./data/venosta_hy1920_swe.nc',
'./data/venosta_hy2021_swe.nc',
'./data/venosta_hy2122_swe.nc',
"./data/venosta_hy1718_swe.nc",
"./data/venosta_hy1819_swe.nc",
"./data/venosta_hy1920_swe.nc",
"./data/venosta_hy2021_swe.nc",
"./data/venosta_hy2122_swe.nc",
]
cog_dirs = [
'./output_cog/venosta_1819_cog/',
'./output_cog/venosta_1819_cog/',
'./output_cog/venosta_1920_cog/',
'./output_cog/venosta_2021_cog/',
'./output_cog/venosta_2122_cog/',
"./output_cog/venosta_1819_cog/",
"./output_cog/venosta_1819_cog/",
"./output_cog/venosta_1920_cog/",
"./output_cog/venosta_2021_cog/",
"./output_cog/venosta_2122_cog/",
]
process_multiple_datasets(net_cdf_paths, cog_dirs)
```
%% Cell type:code id: tags:
``` python
dirs = ["./output_cog/venosta_1718_cog",
"./output_cog/venosta_1819_cog",
"./output_cog/venosta_1920_cog",
"./output_cog/venosta_2021_cog",
"./output_cog/venosta_2122_cog"
]
transform_cogs_in_dirs(dirs, 'EPSG:3035')
# Transform the coordinates of the COGs to any EPSG
# In this case EPSG:3035
dirs = [
"./output_cog/venosta_1718_cog",
"./output_cog/venosta_1819_cog",
"./output_cog/venosta_1920_cog",
"./output_cog/venosta_2021_cog",
"./output_cog/venosta_2122_cog",
]
transform_cogs_in_dirs(dirs, "EPSG:3035")
```
......
swe_maps_stac/transformation.py
+ 18
− 12
View file @ b875c9e3
......@@ -3,6 +3,7 @@ import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from typing import List
def transform_cog(cog_path: str, target_crs: str) -> None:
"""
Transform the coordinates of a COG file.
......@@ -16,16 +17,19 @@ def transform_cog(cog_path: str, target_crs: str) -> None:
"""
with rasterio.open(cog_path) as src:
transform, width, height = calculate_default_transform(
src.crs, target_crs, src.width, src.height, *src.bounds)
src.crs, target_crs, src.width, src.height, *src.bounds
)
kwargs = src.meta.copy()
kwargs.update({
'crs': target_crs,
'transform': transform,
'width': width,
'height': height
})
with rasterio.open('temp.tif', 'w', **kwargs) as dst:
kwargs.update(
{
"crs": target_crs,
"transform": transform,
"width": width,
"height": height,
}
)
with rasterio.open("temp.tif", "w", **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
......@@ -34,9 +38,11 @@ def transform_cog(cog_path: str, target_crs: str) -> None:
src_crs=src.crs,
dst_transform=transform,
dst_crs=target_crs,
resampling=Resampling.nearest)
resampling=Resampling.nearest,
)
os.rename("temp.tif", cog_path)
os.rename('temp.tif', cog_path)
def transform_cogs_in_dirs(dirs: List[str], target_crs: str) -> None:
"""
......@@ -51,6 +57,6 @@ def transform_cogs_in_dirs(dirs: List[str], target_crs: str) -> None:
"""
for dir_path in dirs:
for cog_file in os.listdir(dir_path):
if cog_file.endswith('.tif'):
if cog_file.endswith(".tif"):
cog_path = os.path.join(dir_path, cog_file)
transform_cog(cog_path, target_crs)
swe_maps_stac/utils.py
+ 0
− 1
View file @ b875c9e3
......@@ -41,4 +41,3 @@ def get_cog_coords(ds, row_index, col_index):
"""
x_coord, y_coord = ds.transform * (col_index, row_index)
return x_coord, y_coord
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment