Imagen de geotiff cortada en Python GDAL con archivo geojson

Indagamos en el mundo on line y de este modo regalarte la respuesta para tu inquietud, si continúas con alguna duda puedes dejarnos tu inquietud y te respondemos con mucho gusto.

Solución:

Ahora hay módulos de Python más fáciles de usar para eso, como rasterio

Rasterio emplea GDAL para leer y escribir archivos usando GeoTIFF y muchos otros formatos. Su API utiliza interfaces familiares de Python y SciPy y modismos como administradores de contexto, iteradores y ndarrays.

Por lo tanto de Enmascaramiento de ráster con una entidad poligonal en Recetario Rasterio

import rasterio
from rasterio.mask import mask
# 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("test.tif") as src:
    out_image, out_transform = mask(src, geoms, crop=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)

with rasterio.open("masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

Resultado

Alternativas para archivos

1) Simplemente puede usar los módulos json o geojson para importar la geometría

with open('poly.json') as data_file:    
    geoms= json.load(data_file)
print geoms
u'type': u'Polygon', u'coordinates': [[[250204.0, 141868.0], [250942.0, 141868.0], [250942.0, 141208.0], [250204.0, 141208.0], [250204.0, 141868.0]]]

2) puedes usar el módulo ogr con un shapefile

from osgeo import ogr
reader = ogr.Open('poly.json')
layer = reader.GetLayer()
feature = layer.GetFeature(0)
geoms =json.loads(feature.ExportToJson())['geometry']
print geoms
u'type': u'Polygon', u'coordinates': [[[250204.0, 141868.0], [250942.0, 141868.0], [250942.0, 141208.0], [250204.0, 141208.0], [250204.0, 141868.0]]]

3) también puedes usar el módulo Fiona

Es un contenedor de Python para funciones de acceso a datos vectoriales de la biblioteca OGR

import fiona
with fiona.open("poly.shp") as shapefile:
    geoms = [feature["geometry"] for feature in shapefile]
print geoms
['type': 'Polygon', 'coordinates': [[(250204.0, 141868.0), (250942.0, 141868.0), (250942.0, 141208.0), (250204.0, 141208.0), (250204.0, 141868.0)]]]

Nuevo

Puede usar un filtro como en el script de Luke en ¿Cómo configurar un filtro espacial con Python / GDAL ?. En lugar de cortar, filtra la entrada.

from osgeo import gdal
xmin,ymin,xmax,ymax = [250204.0, 141208.0, 250942.0, 141868.0]
def map_to_pixel(mx,my,gt):
   #Convert from map to pixel coordinates.
   #Only works for geotransforms with no rotation.
   px = int((mx - gt[0]) / gt[1]) #x pixel
   py = int((my - gt[3]) / gt[5]) #y pixel
   return px, py

def extent_to_offset(xmin,ymin,xmax,ymax,gt):
  pxmin,pymin = map_to_pixel(xmin,ymin,gt)
  pxmax,pymax = map_to_pixel(xmax,ymax,gt)
  return pxmin,pymin,abs(pxmax-pxmin),abs(pymax-pymin)

ds=gdal.Open('test.tif')
gt = ds.GetGeoTransform()

bands = ds.RasterCount
band_list = []
#
# Read in bands and store all the data in bandList
#
for i in range(bands):
   band = ds.GetRasterBand(i+1) # 1-based index
   # apply filter to only read the data in the bounding box (xmin, ...)
   data = band.ReadAsArray(xoff, yoff, xsize, ysize)
   band_list.append(data)

 driver = gdal.GetDriverByName("GTiff")
 dst_dst = driver.Create('result.tif', xsize, ysize, 4, gdal.GDT_Byte)
 for j in range(bands):
    data = band_list[j]
    dst_dst.GetRasterBand(j+1).WriteArray(data)

 ....
 dst_dst = None

Solo necesitas agregar las crs

Aquí está mi propia solución. ¡Funciona para cualquier número de bandas, cualquier tipo de geometría (por ejemplo, multipolígono) y funciona con imágenes en cualquier zona!

import geojson as gj
from osgeo import ogr, osr, gdal

# Enable GDAL/OGR exceptions
gdal.UseExceptions()


# GDAL & OGR memory drivers
GDAL_MEMORY_DRIVER = gdal.GetDriverByName('MEM')
OGR_MEMORY_DRIVER = ogr.GetDriverByName('Memory')


def cut_by_geojson(input_file, output_file, shape_geojson):

    # Get coords for bounding box
    x, y = zip(*gj.utils.coords(gj.loads(shape_geojson)))
    min_x, max_x, min_y, max_y = min(x), max(x), min(y), max(y)

    # Open original data as read only
    dataset = gdal.Open(input_file, gdal.GA_ReadOnly)

    bands = dataset.RasterCount

    # Getting georeference info
    transform = dataset.GetGeoTransform()
    projection = dataset.GetProjection()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = -transform[5]

    # Getting spatial reference of input raster
    srs = osr.SpatialReference()
    srs.ImportFromWkt(projection)

    # WGS84 projection reference
    OSR_WGS84_REF = osr.SpatialReference()
    OSR_WGS84_REF.ImportFromEPSG(4326)

    # OSR transformation
    wgs84_to_image_trasformation = osr.CoordinateTransformation(OSR_WGS84_REF,
                                                                srs)
    XYmin = wgs84_to_image_trasformation.TransformPoint(min_x, max_y)
    XYmax = wgs84_to_image_trasformation.TransformPoint(max_x, min_y)

    # Computing Point1(i1,j1), Point2(i2,j2)
    i1 = int((XYmin[0] - xOrigin) / pixelWidth)
    j1 = int((yOrigin - XYmin[1]) / pixelHeight)
    i2 = int((XYmax[0] - xOrigin) / pixelWidth)
    j2 = int((yOrigin - XYmax[1]) / pixelHeight)
    new_cols = i2 - i1 + 1
    new_rows = j2 - j1 + 1

    # New upper-left X,Y values
    new_x = xOrigin + i1 * pixelWidth
    new_y = yOrigin - j1 * pixelHeight
    new_transform = (new_x, transform[1], transform[2], new_y, transform[4],
                     transform[5])

    wkt_geom = ogr.CreateGeometryFromJson(str(shape_geojson))
    wkt_geom.Transform(wgs84_to_image_trasformation)

    target_ds = GDAL_MEMORY_DRIVER.Create('', new_cols, new_rows, 1,
                                          gdal.GDT_Byte)
    target_ds.SetGeoTransform(new_transform)
    target_ds.SetProjection(projection)

    # Create a memory layer to rasterize from.
    ogr_dataset = OGR_MEMORY_DRIVER.CreateDataSource('shapemask')
    ogr_layer = ogr_dataset.CreateLayer('shapemask', srs=srs)
    ogr_feature = ogr.Feature(ogr_layer.GetLayerDefn())
    ogr_feature.SetGeometryDirectly(ogr.Geometry(wkt=wkt_geom.ExportToWkt()))
    ogr_layer.CreateFeature(ogr_feature)

    gdal.RasterizeLayer(target_ds, [1], ogr_layer, burn_values=[1],
                        options=["ALL_TOUCHED=TRUE"])

    # Create output file
    driver = gdal.GetDriverByName('GTiff')
    outds = driver.Create(output_file, new_cols, new_rows, bands,
                          gdal.GDT_Float32)

    # Read in bands and store all the data in bandList
    mask_array = target_ds.GetRasterBand(1).ReadAsArray()
    band_list = []

    for i in range(bands):
        band_list.append(dataset.GetRasterBand(i + 1).ReadAsArray(i1, j1,
                         new_cols, new_rows))

    for j in range(bands):
        data = np.where(mask_array == 1, band_list[j], mask_array)
        outds.GetRasterBand(j + 1).SetNoDataValue(0)
        outds.GetRasterBand(j + 1).WriteArray(data)

    outds.SetProjection(projection)
    outds.SetGeoTransform(new_transform)

    target_ds = None
    dataset = None
    outds = None
    ogr_dataset = None

Aquí puedes ver las comentarios y valoraciones de los lectores

Acuérdate de que te brindamos la opción de interpretar si te fue de ayuda.

---

---

Utiliza Nuestro Buscador