Reputation: 352
I want to:
1
and pixels outside the fields are set to 0
thenCde:
import fiona
import rasterio
import rasterio.mask
import pycrs
def masked_raster(input_file, raster_file):
# Create a masked version of the input raster where pixels falling within one of the fields are set to `1` and pixels outside the fields are set to `0`
data = rasterio.open(raster_file)
#creating the a bounding box with Shapely
## WGS84 coordinates
minx, miny = 24.60, 60.00
maxx, maxy = 25.22, 60.35
bbox = box(minx, miny, maxx, maxy)
#inserting the bounding box into GeoDataFrame
geo = geopandas.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(4326))
#Re-project into the same coordinate system as the raster data
geo = geo.to_crs(crs=data.crs.data)
#get the coordinates of the geometry in a proper format of rasterio
def getFeatures(gdf):
"""Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
import json
return [json.loads(gdf.to_json())['features'][0]['geometry']]
#Getting the geometry coordinates by using the function
coordinates = getFeatures(geo)
# mask to clip the raster with the polygon using the coords variable
out_img, out_transform = mask(data, shapes=coordinates, crop=True)
out_img = rasterio.open(raster_file).read()
return out_img
def reproject_raster(raster_file, dst_crs):
# Reproject the input raster to the provided CRS
src = rasterio.open(raster_file)
# Parse EPSG code
epsg_code = int(data.crs.data['init'][5:])
#copy metadata
out_meta = raster_file.meta.copy()
#update the metadata with new dimensions, transform (affine) and CRS
out_meta.update({"driver": "GTiff", "height": out_img.shape[1], "width": out_img.shape[2], "transform": out_transform, "crs": pycrs.parser.from_epsg_code(epsg_code).to_proj4()})
#save the clipped raster to disk
with rasterio.open(out_tif, "w", **out_meta) as dst:
dst.write(out_img)
dst = src
return dst
# Run this to validate your function works correctly
assert masked_raster('crops.geojson', 'crops.tif')[0].sum() == 1144636.0, "Sorry wrong answer"
assert str(reproject_raster('crops.tif', 'EPSG:4326').crs) == 'EPSG:4326', "Sorry wrong answer"
print("Congratulations, all is working just fine !!!")
Upvotes: 0
Views: 5583
Reputation: 352
Thanks guys, I managed to find the solution to my error: "rasterio.tools.mask.mask (in more recent versions, it is rasterio.mask.mask) includes an option invert. When invert=True, the mask will be applied to pixels that overlap your shape, rather than areas outside the shape"
So I changed:
out_img, out_transform = mask(data, shapes=coordinates, crop=True)
To:
out_img, out_transform = mask(src, geoms, invert=True)
The error is fixed. Here is the full code I rewrote. Though I am not getting the expected output by there is no error, its only the part where it must meet the assertion requirements that must be fixed.
def masked_raster(input_file, raster_file):
# Create a masked version of the input raster where pixels falling within one of the fields are set to `1` and pixels outside the fields are set to `0`
# the polygon GeoJSON geometry
geoms = [{'type': 'Polygon', 'coordinates': [[(250204.0, 141868.0), (250942.0, 141868.0), (250942.0, 141208.0), (250204.0, 141208.0), (250204.0, 141868.0)]]}]
# load the raster, mask it by the polygon and crop it
with rasterio.open("crops.tif") as src:
out_img, out_transform = mask(src, geoms, invert=True)
out_meta = src.meta.copy()
# save the resulting raster
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
return out_img
def reproject_raster(raster_file, dst_crs):
# Reproject the input raster to the provided CRS
with rasterio.open("masked.tif", "w", **out_meta) as dst:
dst.write(out_image)
dst = src
return dst
Upvotes: 5