In [42]:

    

# Import modules
import numpy as np
import gdal
from osgeo import gdal, ogr, gdalnumeric, osr
from PIL import Image
import matplotlib.pyplot as plt
%matplotlib inline
import warnings 
warnings.filterwarnings('ignore')
import os, sys
gdal.UseExceptions()


    

In [43]:

    

# Define raster import function
def import_raster(filename):
    #Open file
    f = gdal.Open(filename)
    md = {'cols':f.RasterXSize}
    md.update({'rows':f.RasterYSize})
    md.update({'n_bands':f.RasterCount})
    md.update({'proj':f.GetGeoTransform()})
    md.update({'xMin':md['proj'][0]})
    md.update({'yMax':md['proj'][3]})
    md.update({'xMax':md['xMin'] + md['cols']/md['proj'][1]})
    md.update({'yMin':md['yMax'] + md['rows']/md['proj'][5]})
    md.update({'extent':(md['xMin'],md['xMax'],md['yMin'],md['yMax'])})
    md.update({'pixelWidth':md['proj'][1]})
    md.update({'pixelHeight':md['proj'][5]})
    
    raster = f.GetRasterBand(1)
    md.update({'noDataVal':raster.GetNoDataValue()})
    md.update({'scaleFactor':raster.GetScale()})

    return(raster,md)


    

In [44]:

    

# Import initial rasters
lidar_f = '/home/annie/Documents/NEON/data/RSDI-2017/NEONdata/SJER/NEON_D17_SJER_DP3_256000_4106000_CHM.tif'
hs_f = '/home/annie/Documents/NEON/data/RSDI-2017/NEONdata/SJER/NEON_D17_SJER_DP3_256000_4106000_CHM.tif'

lidar_rast,lidar_md = import_raster(lidar_f)
hs_rast, hs_md = import_raster(hs_f)


    

In [45]:

    

# Figure out smaller raster
lidar_size = hs_md['rows']*lidar_md['cols']
hs_size = hs_md['rows']*lidar_md['cols']

# this will be the same for the entire area, so grab smaller


    

In [46]:

    

# this allows GDAL to throw Python Exceptions
gdal.UseExceptions()

#  get raster datasource
#
src_ds = gdal.Open(lidar_f)
srcband = src_ds.GetRasterBand(1)

#  create output datasource
dst_layername = "clipping_layer"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )


    

    
Out[46]:





0

In [61]:

    

# Clip raster to created shapefile

# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL.

def clip_raster(shapefile_path, raster_path):
    # Load the source data as a gdalnumeric array
    srcArray = gdalnumeric.LoadFile(raster_path)

    # Also load as a gdal image to get geotransform
    # (world file) info
    srcImage = gdal.Open(raster_path)
    geoTrans = srcImage.GetGeoTransform()

    # Create an OGR layer from a boundary shapefile
    shapef = ogr.Open(shapefile_path)
    lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
    poly = lyr.GetNextFeature()
    print(lyr,poly)

    # Convert the layer extent to image pixel coordinates
    minX, maxX, minY, maxY = lidar_md['extent']
    print(lidar_md['extent'])
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)
    print(ulX,ulY)

    # Calculate the pixel size of the new image
    pxWidth = 1
    pxHeight = 1

    clip = srcArray[:, ulY:lrY, ulX:lrX]

    #
    # EDIT: create pixel offset to pass to new image Projection info
    #
    xoffset =  ulX
    yoffset =  ulY
    print("Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ))

    # Create a new geomatrix for the image
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    # Map points to pixels for drawing the
    # boundary on a blank 8-bit,
    # black and white, mask image.
    points = []
    pixels = []
    geom = poly.GetGeometryRef()
    pts = geom.GetGeometryRef(0)
    for p in range(pts.GetPointCount()):
      points.append((pts.GetX(p), pts.GetY(p)))
    for p in points:
      pixels.append(world2Pixel(geoTrans, p[0], p[1]))
    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
    rasterize = ImageDraw.Draw(rasterPoly)
    rasterize.polygon(pixels, 0)
    mask = imageToArray(rasterPoly)

    # Clip the image using the mask
    clip = gdalnumeric.choose(mask, \
        (clip, 0)).astype(gdalnumeric.uint8)

    # This image has 3 bands so we stretch each one to make them
    # visually brighter
    for i in range(3):
      clip[i,:,:] = stretch(clip[i,:,:])

    # Save new tiff
    #
    #  EDIT: instead of SaveArray, let's break all the
    #  SaveArray steps out more explicity so
    #  we can overwrite the offset of the destination
    #  raster
    #
    ### the old way using SaveArray
    #
    # gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path)
    #
    ###
    #
    gtiffDriver = gdal.GetDriverByName( 'GTiff' )
    if gtiffDriver is None:
        raise ValueError("Can't find GeoTiff Driver")
    gtiffDriver.CreateCopy( "OUTPUT.tif",
        OpenArray( clip, prototype_ds=raster_path, xoff=xoffset, yoff=yoffset )
    )

    gdal.ErrorReset()


    

In [62]:

    

test = clip_raster('/home/annie/Documents/NEON/Day5/group_project/clipping_layer.shp',\
                  hs_f)


    

    




<osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at 0x7f684185f780> > None
(256000.0, 257000.0, 4106000.0, 4107000.0)
0 0






    




---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-62-45b89251fffd> in <module>()
----> 1 test = clip_raster('/home/annie/Documents/NEON/Day5/group_project/clipping_layer.shp',                  hs_f)

<ipython-input-61-0bd44ae34023> in clip_raster(shapefile_path, raster_path)
     30     pxHeight = int(lrY - ulY)
     31 
---> 32     clip = srcArray[:, ulY:lrY, ulX:lrX]
     33 
     34     #

IndexError: too many indices for array

In [64]:

    

world2Pixel?


    

In [ ]: