I am trying to clip a raster using a geoJSON without writing to my disk.
The following code compiles and produces the correct result.
import tempfile
import fiona
import rasterio
from rasterio.io import MemoryFile
import rasterio.mask
from scipy import stats
def getSummStats(raster, shapefile):
#Open Polygon
with fiona.open(shapefile, "r") as shapefile:
shapes = [feature["geometry"] for feature in shapefile]
#Clip the raster
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
out_meta = src.meta
#Setup to save and write the clipped raster in memory
with MemoryFile() as memfile:
with memfile.open(driver = "GTiff", height = out_image.shape[1], width = out_image.shape[2], count=1,
dtype=rasterio.float64, transform = out_transform) as dataset:
#write the clipped raster to memory
dataset.write(out_image)
#read the clipped raster from memory
ds = dataset.read()
return stats.describe(ds, axis = None)
with rasterio.open("B:/dd.weather.gc.ca/model_gem_global/15km/grib2/lat_lon/00/000/CMC_glb_TMP_TGL_2_latlon.15x.15_2020052500_P000_Copy.grib2") as src:
print(getSummStats(src, "C:/Users/Tom/Downloads/canada_provinces/canada_provinces/canada_provinces.geojson"))
However, I would like to be able to call this script from the command line and input the raster and geoJSON when the script is called. I made some changes to accomplish this.
import tempfile
import fiona
import rasterio
from rasterio.io import MemoryFile
import rasterio.mask
from scipy import stats
import os, sys
def getSummStats(raster, geoJSON_path):
#Open Polygon
with fiona.open(geoJSON_path, "r") as shapefile:
shapes = [feature["geometry"] for feature in shapefile]
#Clip the raster
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
out_meta = src.meta
#Setup to save and write the clipped raster in memory
with MemoryFile() as memfile:
with memfile.open(driver = "GTiff", height = out_image.shape[1], width = out_image.shape[2], count=1,
dtype=rasterio.float64, transform = out_transform) as dataset:
#write the clipped raster to memory
dataset.write(out_image)
#read the clipped raster from memory
ds = dataset.read()
return stats.describe(ds, axis = None)
def main(raster_path, geoJSON_path):
with rasterio.open(raster_path) as src:
print(getSummStats(src, geoJSON_path))
if name == 'main':
#
# example run : $ Polygon_CMDVersion.py /<full-path>/<raster-name>.grib2 /<full-path>/<geoJSON-name>.geoJSON
#
if len( sys.argv ) < 2:
print ("[ ERROR ] you must two args. 1) the full raster path and 2) the full geoJSON path")
sys.exit( 1 )
main(sys.argv[1], sys.argv[2])
However, this second method gives me the following error
File "C:\Users\Tom\Downloads\Polygon_CMDVersion.py" line 21, in getSummStats
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
NameError: name 'src' is not defined
If this helps, here is what I wrote in the command line
Polygon_CMDVersion.py B:/dd.weather.gc.ca/model_gem_global/15km/grib2/lat_lon/00/000/CMC_glb_TMP_TGL_2_latlon.15x.15_2020052500_P000_Copy.grib2 C:/Users/Tom/Downloads/canada_provinces/canada_provinces/canada_provinces.geoJSON
How can I resolve this issue?
def getSummStats(raster, geoJSON_path):todef getSummStats(src, geoJSON_path):? – Johan Jun 15 '20 at 15:09srcin yourgetSummStatsfunction, nor is it passed in as an argument. You have arasterarg but do not use it. – user2856 Jun 15 '20 at 23:38