gdalpythonqgis-3rasterioshapefile
I clipped my TIFF image to smaller TIFF via using a shapefile as mask with rasterio Python library. After that my square shaped clip has bigger extent than shapefile and outside the the shapefile area is black and valued as 0. How can I change it to NODATA values?
I took one png file from the data_dir of GeoServer. Install GeoServer and you will have the same image to play with in directory:
\geoserver-2.7.1\data_dir\coverages\mosaic_sample\
Convert png into tiff and assign projection.
gdal_translate -a_srs epsg:4326 global_mosaic_6.png tiff.tiff
Input file size is 50, 50
0...10...20...30...40...50...60...70...80...90...100 - done.
Then use gdalwarp:
gdalwarp -cutline burn.shp -crop_to_cutline tiff.tiff cut.tif
Creating output file that is 16P x 15L.
Processing input file tiff.tiff.
Warning : the source raster dataset has a SRS, but the cutline features
not. We assume that the cutline coordinates are expressed in the destination SRS. If not, cutline results may be incorrect.
0...10...20...30...40...50...60...70...80...90...100 - done.
This is the result:
Burn.shp contains one polygon:
POLYGON (( 10.375610236185853 39.384332051098106, 10.344426149120759 39.730475417520644, 10.87455562922735 39.92069834861771, 11.152094004106685 39.546489303836594, 10.92133175982499 39.16292503293594, 10.375610236185853 39.384332051098106 ))
It appeared later from the comments that the cutline shapefile "polygon border200millas_polygon.shp" contained three polygons and gdalwarp seems to accept only one geometry as a cutline. However, it can take a multipolygon which can be make from polygons for example with the QGIS function "dissolve".
The issue is resolved. The issue was I misread the documentation. On a second read, the rasterio.mask documentation clearly states that polygons should be a list of GeoJSON-like dicts. I found the following snippet of code to generate those lists from this answer:
geoms = [feature["geometry"] for feature in shapefile]
Here is the the full code that is working:
import rasterio
from rasterio import features
from rasterio import mask
import json
import geopandas as gpd
import os
results = []
final_results = []
with rasterio.open(r"C:\1927_oahu\tif\_Line1_6to8_0.tif") as src:
src_meta = src.meta.copy()
src_affine = src_meta.get("transform")
band = src.read(1)
for geometry, raster_value in features.shapes(band, transform=src_affine):
if raster_value == 1:
result = {'properties': {'raster_value': raster_value}, 'geometry': geometry}
results.append(result)
gpd_results = gpd.GeoDataFrame.from_features(results)
gpd_results["area"] = gpd_results["geometry"].area
gpd_results_filtered = gpd_results[gpd_results["area"] > 5000]
gpd_filtered_json_str = gpd_results_filtered.to_json()
gpd_filtered_json_dict = json.loads(gpd_filtered_json_str)
final_results = [feature["geometry"] for feature in gpd_filtered_json_dict["features"]]
masked_band, masked_transform = mask.mask(src, final_results, invert=True)
masked_band[masked_band > 255] = 0
src_meta.update(dtype=rasterio.uint8, height=int(masked_band.shape[1]), width=int(masked_band.shape[2]), nodata=0, transform=masked_transform, compress='lzw')
with rasterio.open(r"C:\1927_oahu\output\_Line1_6to8_0.tif"), 'w', **src_meta) as dst:
dst.write(masked_band.astype(rasterio.uint8))
Best Answer
You need to open the source file to read data values and metadata, then read the band data as numpy array, make the required changes to the numpy array and then save numpay array to a new dataset file:
References:
https://rasterio.groups.io/g/main/topic/change_the_nodata_value_in_a/28801885 https://rasterio.readthedocs.io/en/latest/topics/writing.html