0

I want to:

  1. 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 then
  2. Reproject the input raster to the provided CRS

Cde:

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 !!!")

Error: enter image description here

Git Repo Link

1
  • Please proofread after posting (edited). And what exactly has this to do with machine-learning? Commented Mar 14, 2021 at 19:07

1 Answer 1

5

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
Sign up to request clarification or add additional context in comments.

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.