Guest User

Untitled

a guest
Jan 18th, 2019
89
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.46 KB | None | 0 0
  1. with fiona.open("C:/.../poly.shp", "r") as shapefile:
  2. shapes = [feature["geometry"] for feature in shapefile]
  3.  
  4. with rasterio.open("D:/.../MODIS_EVI_2006.tif") as src:
  5. out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
  6. out_meta = src.meta.copy()
  7.  
  8. out_meta.update({"driver": "GTiff",
  9. "height": out_image.shape[1],
  10. "width": out_image.shape[2],
  11. "transform": out_transform})
  12.  
  13. with rasterio.open("mask2006.tif", "w", **out_meta) as dest:
  14. dest.write(out_image)
  15.  
  16. from osgeo import gdal
  17. from osgeo import ogr
  18. from osgeo import gdalconst
  19.  
  20. ndsm = 'D:/../MODIS_EVI_2000.tif'
  21. shp = 'C:/Doy2000/Pivos/1996poly.shp'
  22. data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
  23. geo_transform = data.GetGeoTransform()
  24. #source_layer = data.GetLayer()
  25. x_min = geo_transform[0]
  26. y_max = geo_transform[3]
  27. x_max = x_min + geo_transform[1] * data.RasterXSize
  28. y_min = y_max + geo_transform[3] * data.RasterYSize
  29. x_res = data.RasterXSize
  30. y_res = data.RasterYSize
  31. mb_v = ogr.Open(shp)
  32. mb_l = mb_v.GetLayer()
  33. pixel_width = geo_transform[1]
  34. output = 'C:/Doy2000/Pivos/my2000.tif'
  35. target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
  36. target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
  37. band = target_ds.GetRasterBand(1)
  38. NoData_value = -999999
  39. band.SetNoDataValue(NoData_value)
  40. band.FlushCache()
  41. gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow"])
  42.  
  43. target_ds = None
Add Comment
Please, Sign In to add comment