Geographic Information Systems Asked by FelixB. on February 9, 2021
I want to subset multiple raster files with multiple shapefiles and the following code is working (it is not the whole code, only the part for the subsetting).
I am using Python 3.6.8
i_shp = 0 # iterate over .shp-files
while i_shp < len(shp_list):
with fiona.open(shp_list[i_shp], "r") as shapefile:
shapes = [feature["geometry"] for feature in shapefile]
i_ras = 0 # iterate over raster files
while i_ras < len(raster_list):
for scene in raster_list:
with rasterio.open(raster_list[i_ras], "r") as src:
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
out_meta = src.meta
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
with rasterio.open(
str(outpath) + raster_names[i_ras] + str(shp_names[i_shp]) + str(ras_extension[1:])
, "w", **out_meta) as dest:
dest.write(out_image)
i_ras = i_ras + 1
i_shp = i_shp + 1
I`ve got 7 rasters and 2 shapefiles and after I run the code I have 14 subsets.
But I want to place this code inside a function but then it’s not working correctly. I just get 7 subsets (only from the first shapefile in my shp_list
).
I think that the while-loop isn’t working correctly, but I cannot figure out why.
def subs(shp_list, shp_names, raster_list, raster_names, outpath):
# subsetting all raster files with all .shp-files
i_shp = 0 # iterate over .shp-files
while i_shp < len(shp_list):
with fiona.open(shp_list[i_shp], "r") as shapefile:
shapes = [feature["geometry"] for feature in shapefile]
i_ras = 0 # iterate over raster files
while i_ras < len(raster_list):
for scene in raster_list:
with rasterio.open(raster_list[i_ras], "r") as src:
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
out_meta = src.meta
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
with rasterio.open(
str(outpath) + raster_names[i_ras] + str(shp_names[i_shp]) + str(ras_extension[1:])
, "w", **out_meta) as dest:
dest.write(out_image)
i_ras = i_ras + 1
i_shp = i_shp + 1
# number of created subsets
subset_count = i_shp * i_ras
if len(shp_list) * len(raster_list) == subset_count:
(str("Done. n ") + str(subset_count) + str(" subsets created"))
return print(str("Done. n ") + str(subset_count) + str(" subsets created"))
You can use the script below. You don't need while
here. Instead, use for .... in enumerate(...):
for a cleaner script. Also I've changed string definitions to f-strings
which work in Python 3.6 and provide a clenaer string definition.
# other imports
from rasterio.mask import mask
# other stuffs -> shp_list, shp_names etc.
def subs(shp_list, shp_names, raster_list, raster_names, outpath):
for i, shp in enumerate(shp_list):
with fiona.open(shp, "r") as shapefile:
shapes = [feature["geometry"] for feature in shapefile]
for j, ras in enumerate(raster_list):
for scene in raster_list:
with rasterio.open(ras, "r") as src:
out_image, out_transform = mask(src, shapes, crop=True)
out_meta = src.meta
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
ras_path = f"{outpath}{raster_names[j]}{shp_names[i]}{ras_extension[1:]}"
with rasterio.open(ras_path, "w", **out_meta) as dest:
dest.write(out_image)
# number of created subsets
subset_count = (i+1) * (j+1)
if len(shp_list) * len(raster_list) == subset_count:
print(f"Done. n{subset_count} subsets created")
# calling the function
subs(shp_list, shp_names, raster_list, raster_names, outpath)
Answered by Kadir Şahbaz on February 9, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP