Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Crs checking #9

Merged
merged 2 commits into from
Feb 28, 2024
Merged

Crs checking #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a49d9bd
improve crs handling
mayrajeo Feb 28, 2024
5133272
add library
mayrajeo Feb 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 10 additions & 6 deletions geo2ml/data/tabular.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,9 +56,11 @@ def sample_raster_with_points(
gdf = gpd.read_file(sampling_locations, layer=gpkg_layer)

with rio.open(input_raster) as src:
assert str(gdf.crs) == str(
src.crs
), "Sampling locations and input raster have different crs"
if src.gcps[1]:
in_crs = src.gcps[1]
else:
in_crs = src.crs
gdf = gdf.to_crs(in_crs)
coords = [(x, y) for x, y in zip(gdf.geometry.x, gdf.geometry.y)]
values = np.array([p for p in src.sample(coords)])
prof = src.profile
Expand Down Expand Up @@ -104,9 +106,11 @@ def sample_raster_with_polygons(

gdf = gpd.read_file(sampling_locations, layer=gpkg_layer)
with rio.open(input_raster) as src:
assert str(gdf.crs) == str(
src.crs
), "Sampling locations and input raster have different crs"
if src.gcps[1]:
in_crs = src.gcps[1]
else:
in_crs = src.crs
gdf = gdf.to_crs(in_crs)
n_bands = src.count
prof = src.profile
zstats = []
Expand Down
22 changes: 15 additions & 7 deletions geo2ml/data/tiling.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,13 +86,17 @@ def tile_raster(
(dy, dy + self.gridsize_y), (dx, dx + self.gridsize_x)
)
prof = src.profile.copy()

if src.gcps[1]:
in_crs = src.gcps[1]
else:
in_crs = src.crs
prof.update(
height=window.height,
width=window.width,
transform=rio_windows.transform(window, src.transform),
compress="lzw",
predictor=2,
crs=in_crs,
)
fname = f"R{iy}C{ix}"
data = src.read(window=window)
Expand All @@ -104,9 +108,8 @@ def tile_raster(
dest.write(data)
polys.append(box(*dest.bounds))
names.append(fname)
self.grid = gpd.GeoDataFrame(
{"cell": names, "geometry": polys}, crs=src.crs
)

self.grid = gpd.GeoDataFrame({"cell": names, "geometry": polys}, crs=in_crs)
return

def tile_vector(
Expand Down Expand Up @@ -139,7 +142,7 @@ def tile_vector(
raise Exception("Unknown output format, must be either `geojson` or `gpkg`")

vector = gpd.read_file(path_to_vector, layer=gpkg_layer)

vector = vector.to_crs(self.grid.crs)
sindex = vector.sindex
for row in tqdm(self.grid.itertuples()):
possible_matches_index = list(sindex.intersection(row.geometry.bounds))
Expand Down Expand Up @@ -198,14 +201,15 @@ def tile_and_rasterize_vector(
os.makedirs(self.rasterized_vector_path)

vector = gpd.read_file(path_to_vector)
vector = vector.to_crs(self.grid.crs)
le = LabelEncoder()
vector["label"] = (
le.fit_transform(vector[column].values) + 1
) # We want the labels to start from 1 as 0 is background most of the time
with open(self.outpath / "label_map.txt", "w") as f:
for c, i in zip(le.classes_, le.transform(le.classes_)):
f.write(f"{c}: {i+1}\n")
for r in tqdm(os.listdir(self.raster_path)):
for r in tqdm([f for f in os.listdir(self.raster_path) if f.endswith("tif")]):
with rio.open(self.raster_path / r) as src:
out_meta = src.meta
out_meta.update({"count": 1})
Expand Down Expand Up @@ -248,6 +252,10 @@ def untile_raster(

for f in rasters:
src = rio.open(f)
if src.gcps[1]:
in_crs = src.gcps[1]
else:
in_crs = src.crs
files_to_mosaic.append(src)

mosaic, out_tfm = rio_merge(files_to_mosaic, method=method)
Expand All @@ -259,7 +267,7 @@ def untile_raster(
"height": mosaic.shape[1],
"width": mosaic.shape[2],
"transform": out_tfm,
"crs": src.crs,
"crs": in_crs,
}
)

Expand Down
8 changes: 6 additions & 2 deletions nbs/10_data.tabular.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,9 @@
" gdf = gpd.read_file(sampling_locations, layer=gpkg_layer)\n",
"\n",
" with rio.open(input_raster) as src:\n",
" assert str(gdf.crs) == str(src.crs), \"Sampling locations and input raster have different crs\"\n",
" if src.gcps[1]: in_crs = src.gcps[1]\n",
" else: in_crs = src.crs\n",
" gdf = gdf.to_crs(in_crs)\n",
" coords = [(x,y) for x,y in zip(gdf.geometry.x, gdf.geometry.y)]\n",
" values = np.array([p for p in src.sample(coords)])\n",
" prof = src.profile\n",
Expand Down Expand Up @@ -793,7 +795,9 @@
" \n",
" gdf = gpd.read_file(sampling_locations, layer=gpkg_layer)\n",
" with rio.open(input_raster) as src:\n",
" assert str(gdf.crs) == str(src.crs), \"Sampling locations and input raster have different crs\"\n",
" if src.gcps[1]: in_crs = src.gcps[1]\n",
" else: in_crs = src.crs\n",
" gdf = gdf.to_crs(in_crs)\n",
" n_bands = src.count\n",
" prof = src.profile\n",
" zstats = []\n",
Expand Down
48 changes: 27 additions & 21 deletions nbs/12_data.tiling.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,13 +234,15 @@
" window = rio_windows.Window.from_slices((dy, dy+self.gridsize_y),\n",
" (dx, dx+self.gridsize_x))\n",
" prof = src.profile.copy()\n",
" \n",
" if src.gcps[1]: in_crs = src.gcps[1]\n",
" else: in_crs = src.crs\n",
" prof.update(\n",
" height=window.height,\n",
" width=window.width,\n",
" transform= rio_windows.transform(window, src.transform),\n",
" compress='lzw',\n",
" predictor=2\n",
" predictor=2,\n",
" crs=in_crs\n",
" )\n",
" fname = f'R{iy}C{ix}'\n",
" data = src.read(window=window)\n",
Expand All @@ -252,7 +254,8 @@
" dest.write(data)\n",
" polys.append(box(*dest.bounds))\n",
" names.append(fname)\n",
" self.grid = gpd.GeoDataFrame({'cell': names, 'geometry':polys}, crs=src.crs) \n",
"\n",
" self.grid = gpd.GeoDataFrame({'cell': names, 'geometry':polys}, crs=in_crs) \n",
" return\n",
" \n",
" def tile_vector(self, path_to_vector:Path|str, min_area_pct:float=0.0, gpkg_layer:str=None, output_format:str='geojson') -> None:\n",
Expand Down Expand Up @@ -280,7 +283,7 @@
" )\n",
" \n",
" vector = gpd.read_file(path_to_vector, layer=gpkg_layer)\n",
" \n",
" vector = vector.to_crs(self.grid.crs)\n",
" sindex = vector.sindex\n",
" for row in tqdm(self.grid.itertuples()):\n",
" possible_matches_index = list(sindex.intersection(row.geometry.bounds))\n",
Expand Down Expand Up @@ -322,12 +325,13 @@
" if not os.path.exists(self.rasterized_vector_path): os.makedirs(self.rasterized_vector_path)\n",
"\n",
" vector = gpd.read_file(path_to_vector)\n",
" vector = vector.to_crs(self.grid.crs)\n",
" le = LabelEncoder()\n",
" vector['label'] = le.fit_transform(vector[column].values) + 1 # We want the labels to start from 1 as 0 is background most of the time\n",
" with open(self.outpath/'label_map.txt', 'w') as f:\n",
" for c, i in zip(le.classes_, le.transform(le.classes_)):\n",
" f.write(f'{c}: {i+1}\\n')\n",
" for r in tqdm(os.listdir(self.raster_path)):\n",
" for r in tqdm([f for f in os.listdir(self.raster_path) if f.endswith('tif')]):\n",
" with rio.open(self.raster_path/r) as src:\n",
" out_meta = src.meta\n",
" out_meta.update({'count':1})\n",
Expand Down Expand Up @@ -415,7 +419,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "d8375a086cab444694f9820081d1cc51",
"model_id": "51eaf1ee31ff4becaca80ee5acd7a233",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -441,7 +445,7 @@
"text/markdown": [
"---\n",
"\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L111){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L112){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"\n",
"### Tiler.tile_vector\n",
"\n",
Expand All @@ -455,7 +459,7 @@
"text/plain": [
"---\n",
"\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L111){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L112){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"\n",
"### Tiler.tile_vector\n",
"\n",
Expand Down Expand Up @@ -503,7 +507,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "78a8237ecfe2400aa44d4ad5e987a245",
"model_id": "36ab9db52432439ab7226f8c49d249f0",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -527,7 +531,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "1b8e2bbc8a1e4cb798dfdcd5ab14e783",
"model_id": "5eef8bacd5ca427ba2ce5da0b52fd34f",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -562,7 +566,7 @@
"text/markdown": [
"---\n",
"\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L171){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L174){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"\n",
"### Tiler.tile_and_rasterize_vector\n",
"\n",
Expand All @@ -577,7 +581,7 @@
"text/plain": [
"---\n",
"\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L171){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"[source](https://github.com/mayrajeo/geo2ml/blob/main/geo2ml/data/tiling.py#L174){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n",
"\n",
"### Tiler.tile_and_rasterize_vector\n",
"\n",
Expand Down Expand Up @@ -619,7 +623,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "8d84a2c087094c8890b0c0aadc011ad2",
"model_id": "c181d67bfe8a4ff3b5c6661d42b71788",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -644,7 +648,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "08264d32474d430790dec1232bbd4d30",
"model_id": "1620997e6e7049488536406b7181dcad",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -718,6 +722,8 @@
" \n",
" for f in rasters:\n",
" src = rio.open(f)\n",
" if src.gcps[1]: in_crs = src.gcps[1]\n",
" else: in_crs = src.crs\n",
" files_to_mosaic.append(src)\n",
" \n",
" mosaic, out_tfm = rio_merge(files_to_mosaic, method=method)\n",
Expand All @@ -727,7 +733,7 @@
" 'height': mosaic.shape[1],\n",
" 'width': mosaic.shape[2],\n",
" 'transform': out_tfm,\n",
" 'crs': src.crs})\n",
" 'crs': in_crs})\n",
" \n",
" with rio.open(outfile, 'w', **out_meta) as dest: dest.write(mosaic)\n",
" \n",
Expand Down Expand Up @@ -783,7 +789,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "6816c16c29cb483e9f81ab6adb84818e",
"model_id": "28ffa470c0b6451fb0fa4c232bdd2351",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -814,7 +820,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "2ade4ef218d24dfca1ad7f5d48470b2b",
"model_id": "4ce5f76ad1634f2b8d1035398ea14aee",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -893,7 +899,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "aae41955234c4dbb9c343e832f1019be",
"model_id": "9f2f690239464a208a26423622586de6",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -918,7 +924,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "31c511247bd647b29960408c71852eaf",
"model_id": "d309e28e01b64a6eb4f2e91bc1d5765a",
"version_major": 2,
"version_minor": 0
},
Expand All @@ -932,7 +938,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "0901b05ad37040ef81dad72ff9089d42",
"model_id": "3d90ef21fc42458f8185f86312299656",
"version_major": 2,
"version_minor": 0
},
Expand Down Expand Up @@ -965,7 +971,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "cc0696a910f349cea9e468ab7c852bf0",
"model_id": "54a784a96e4b408e8e72760f55a1f252",
"version_major": 2,
"version_minor": 0
},
Expand Down
Binary file modified nbs/example_data/R70C21.dbf
View file
Open in desktop
Binary file not shown.
2 changes: 1 addition & 1 deletion nbs/example_data/tiles/coco_norm.json
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion nbs/example_data/tiles/coco_rot.json
View file
Open in desktop

Large diffs are not rendered by default.

Binary file modified nbs/example_data/tiles/vectors.gpkg
View file
Open in desktop
Binary file not shown.
Binary file modified nbs/example_data/tiles_partial/vectors.gpkg
View file
Open in desktop
Binary file not shown.
Loading