User Tools

Site Tools


wiki:python:basicpython3

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

wiki:python:basicpython3 [2018/06/01 20:13] (current)
Line 1: Line 1:
 +
 +===== Some advanced modules =====
 +
 +<code python | 07/​01-shaded-relief.py>​
 +    """​
 +    Creates a shaded relief ASCII grid
 +    from an ASCII DEM.  Also outputs
 +    intermediate grids for slope and
 +    aspect.
 +    """​
 +    from linecache import getline
 +    import numpy as np
 +     
 +    # File name of ASCII digital elevation model
 +    source = "​dem.asc"​
 +    # File name of the slope grid
 +    slopegrid = "​slope.asc"​
 +    # File name of the aspect grid
 +    aspectgrid = "​aspect.asc"​
 +    # Output file name for shaded relief
 +    shadegrid = "​relief.asc"​
 +    ## Shaded elevation parameters
 +    # Sun direction
 +    azimuth=315.0
 +    # Sun angle
 +    altitude=45.0
 +    # Elevation exageration
 +    z=1.0
 +    # Resolution
 +    scale=1.0
 +    # No data value for output
 +    NODATA = -9999
 +     
 +    # Needed for numpy conversions
 +    deg2rad = 3.141592653589793 / 180.0
 +    rad2deg = 180.0 / 3.141592653589793
 +     
 +    # Parse the header using a loop and
 +    # the built-in linecache module
 +    hdr = [getline(source,​ i) for i in range(1,7)]
 +    values = [float(h.split("​ "​)[-1].strip()) for h in hdr]
 +    cols,​rows,​lx,​ly,​cell,​nd = values
 +    xres = cell
 +    yres = cell * -1
 +     
 +    # Load the dem into a numpy array
 +    arr = np.loadtxt(source,​ skiprows=6)
 +     
 +    # Exclude 2 pixels around the edges which are usually NODATA.
 +    # Also set up structure for a 3x3 windows to process the slope
 +    # throughout the grid
 +    window = []
 +    for row in range(3):
 +        for col in range(3):
 +            window.append(arr[row:​(row + arr.shape[0] - 2), \
 +            col:(col + arr.shape[1] - 2)])
 +     
 +    # Process each cell
 +    x = ((z * window[0] + z * window[3] + z * window[3] + z * window[6]) \
 +       - (z * window[2] + z * window[5] + z * window[5] + z * window[8])) \
 +      / (8.0 * xres * scale);
 +     
 +    y = ((z * window[6] + z * window[7] + z * window[7] + z * window[8]) \
 +       - (z * window[0] + z * window[1] + z * window[1] + z * window[2])) \
 +      / (8.0 * yres * scale);
 +     
 +    # Calculate slope  ​
 +    slope = 90.0 - np.arctan(np.sqrt(x*x + y*y)) * rad2deg
 +     
 +    # Calculate aspect
 +    aspect = np.arctan2(x,​ y)  ​
 +     
 +    # Calculate the shaded relief
 +    shaded = np.sin(altitude * deg2rad) * np.sin(slope * deg2rad) \
 +           + np.cos(altitude * deg2rad) * np.cos(slope * deg2rad) \
 +           * np.cos((azimuth - 90.0) * deg2rad - aspect);
 +    shaded = shaded * 255
 +     
 +    # Rebuild the new header
 +    header = "​ncols ​       %s\n" % shaded.shape[1]
 +    header += "​nrows ​       %s\n" % shaded.shape[0]
 +    header += "​xllcorner ​   %s\n" % (lx + (cell * (cols - shaded.shape[1])))
 +    header += "​yllcorner ​   %s\n" % (ly + (cell * (rows - shaded.shape[0])))
 +    header += "​cellsize ​    ​%s\n"​ % cell
 +    header += "​NODATA_value ​     %s\n" % NODATA
 +     
 +    # Set no-data values
 +    for pane in window:
 +       ​slope[pane == nd] = NODATA
 +       ​aspect[pane == nd] = NODATA
 +       ​shaded[pane == nd] = NODATA
 +     
 +    # Open the output file, add the header, save the slope grid
 +    with open(slopegrid,​ "​wb"​) as f:
 +      f.write(header)
 +      np.savetxt(f,​ slope, fmt="​%4i"​)
 +     
 +    # Open the output file, add the header, save the slope grid
 +    with open(aspectgrid,​ "​wb"​) as f:
 +      f.write(header)
 +      np.savetxt(f,​ aspect, fmt="​%4i"​)
 +     
 +    # Open the output file, add the header, save the array
 +    with open(shadegrid,​ "​wb"​) as f:
 +      f.write(header)
 +      np.savetxt(f,​ shaded, fmt="​%4i"​)
 +
 +</​code>​
 +
 +
 +<code python | 07/​02-contour.py>​
 +    """​
 +    Use GDAL and OGR to create a contour shapefile
 +    """​
 +    import gdal
 +    import ogr
 +     
 +    # Elevation DEM
 +    source = "​dem.asc"​
 +    # Output shapefile
 +    target = "​contour"​
 +     
 +     
 +    ogr_ds = ogr.GetDriverByName('​ESRI Shapefile'​).CreateDataSource(target + "​.shp"​)
 +    ogr_lyr = ogr_ds.CreateLayer(target,​ geom_type = ogr.wkbLineString25D)
 +    field_defn = ogr.FieldDefn('​ID',​ ogr.OFTInteger)
 +    ogr_lyr.CreateField(field_defn)
 +    field_defn = ogr.FieldDefn('​ELEV',​ ogr.OFTReal)
 +    ogr_lyr.CreateField(field_defn)
 +     
 +    # gdal.ContourGenerate() arguments
 +    # Band srcBand,
 +    # double contourInterval,​
 +    # double contourBase,​
 +    # double[] fixedLevelCount,​
 +    # int useNoData,
 +    # double noDataValue,​
 +    # Layer dstLayer,
 +    # int idField,
 +    # int elevField
 +     
 +    ds = gdal.Open('​dem.asc'​)
 +    gdal.ContourGenerate(ds.GetRasterBand(1),​ 400, 10, [], 0, 0, ogr_lyr, 0, 1)
 +
 +
 +</​code>​
 +
 +
 +<code python | 07/​03-drawContours.py>​
 +    """​
 +    Draw an entire contour shapefile
 +    to a pngcanvas image.
 +    """​
 +    import shapefile
 +    import pngcanvas
 +    # Open the contours
 +    r = shapefile.Reader("​contour.shp"​)
 +    # Setup the world to pixels conversion
 +    xdist = r.bbox[2] - r.bbox[0]
 +    ydist = r.bbox[3] - r.bbox[1]
 +    iwidth = 800
 +    iheight = 600
 +    xratio = iwidth/​xdist
 +    yratio = iheight/​ydist
 +    contours = []
 +    # Loop through all shapes
 +    for shape in r.shapes():
 +      # Loop through all parts
 +      for i in range(len(shape.parts)):​
 +        pixels=[]
 +        pt = None
 +        if i<​len(shape.parts)-1:​
 +          pt = shape.points[shape.parts[i]:​shape.parts[i+1]]
 +        else:
 +          pt = shape.points[shape.parts[i]:​]
 +        for x,y in pt:
 +          px = int(iwidth - ((r.bbox[2] - x) * xratio))
 +          py = int((r.bbox[3] - y) * yratio)
 +          pixels.append([px,​py])
 +        contours.append(pixels)
 +    # Set up the output canvas
 +    canvas = pngcanvas.PNGCanvas(iwidth,​iheight)
 +    # PNGCanvas accepts rgba byte arrays for colors
 +    red = [0xff,​0,​0,​0xff]
 +    canvas.color = red 
 +    # Loop through the polygons and draw them
 +    for c in contours:
 +      canvas.polyline(c)
 +    # Save the image
 +    f = open("​contours.png",​ "​wb"​)
 +    f.write(canvas.dump())
 +    f.close()
 +
 +
 +</​code>​
 +
 +
 +
 +
 +<code python | 07/​04-grid.py>​
 +    """​
 +    Converts a LIDAR LAS file to an
 +    ASCII DEM.  Interpolation is used
 +    to account for data loss.  ​
 +    """​
 +    # We only need 2 modules
 +    # both available on PyPI!
 +    from laspy.file import File
 +    import numpy as np
 +     
 +    # Source LAS file
 +    source = "​lidar.las"​
 +     
 +    # Output ASCII DEM file
 +    target = "​lidar.asc"​
 +     
 +    # Grid cell size (data units)
 +    cell = 1.0
 +     
 +    # No data value for output DEM
 +    NODATA = 0
 +     
 +    # Open LIDAR LAS file
 +    las = File(source,​ mode="​r"​)
 +     
 +    #xyz min and max
 +    min = las.header.min
 +    max = las.header.max
 +     
 +    # Get the x axis distance
 +    xdist = max[0] - min[0]
 +     
 +    # Get the y axis distance
 +    ydist = max[1] - min[1]
 +     
 +    # Number of columns for our grid
 +    cols = int(xdist) / cell
 +     
 +    # Number of rows for our grid
 +    rows = int(ydist) / cell
 +     
 +    cols += 1
 +    rows += 1
 +     
 +    # Track how many elevation
 +    # values we aggregate
 +    count = np.zeros((rows,​ cols)).astype(np.float32)
 +    # Aggregate elevation values
 +    zsum = np.zeros((rows,​ cols)).astype(np.float32)
 +     
 +    # Y resolution is negative
 +    ycell = -1 * cell
 +     
 +    # Project x,y values to grid
 +    projx = (las.x - min[0]) / cell
 +    projy = (las.y - min[1]) / ycell
 +    # Cast to integers and clip for use as index
 +    ix = projx.astype(np.int32)
 +    iy = projy.astype(np.int32)
 +     
 +    # Loop through x,y,z arrays, add to grid shape,
 +    # and aggregate values for averaging
 +    for x,y,z in np.nditer([ix,​ iy, las.z]):
 +      count[y, x]+=1
 +      zsum[y, x]+=z
 +     
 +    # Change 0 values to 1 to avoid numpy warnings, ​
 +    # and NaN values in array
 +    nonzero = np.where(count>​0,​ count, 1)
 +    # Average our z values
 +    zavg = zsum/​nonzero
 +     
 +    # Interpolate 0 values in array to avoid any
 +    # holes in the grid
 +    mean = np.ones((rows,​cols)) * np.mean(zavg)
 +    left = np.roll(zavg,​ -1, 1)
 +    lavg = np.where(left>​0,​left,​mean)
 +    right = np.roll(zavg,​ 1, 1)
 +    ravg = np.where(right>​0,​right,​mean)
 +    interpolate = (lavg+ravg)/​2
 +    fill=np.where(zavg>​0,​zavg,​interpolate)
 +     
 +    # Create our ASCII DEM header
 +    header = "​ncols ​       %s\n" % fill.shape[1]
 +    header += "​nrows ​       %s\n" % fill.shape[0]
 +    header += "​xllcorner ​   %s\n" % min[0]
 +    header += "​yllcorner ​   %s\n" % min[1]
 +    header += "​cellsize ​    ​%s\n"​ % cell
 +    header += "​NODATA_value ​     %s\n" % NODATA
 +     
 +    # Open the output file, add the header, save the array
 +    with open(target,​ "​wb"​) as f:
 +      f.write(header)
 +      # The fmt string ensures we output floats ​
 +      # that have at least one number but only 
 +      # two decimal places
 +      np.savetxt(f,​ fill, fmt="​%1.2f"​)
 +
 +</​code>​
 +
 +
 +<code python | 07/​05-dem2img.py>​
 +    """​Convert an ASCII DEM to an image."""​
 +    import numpy as np
 +    import Image
 +    import ImageOps
 +     
 +    # Source LIDAR DEM file
 +    source = "​relief.asc"​
 +     
 +    # Output image file
 +    target = "​relief.bmp"​
 +     
 +    # Load the ASCII DEM into a numpy array
 +    arr = np.loadtxt(source,​ skiprows=6)
 +     
 +    # Convert array to numpy image
 +    im = Image.fromarray(arr).convert('​RGB'​)
 +     
 +    # Enhance the image:
 +    # equalize and increase contrast
 +    im = ImageOps.equalize(im)
 +    im = ImageOps.autocontrast(im)
 +     
 +    # Save the image
 +    im.save(target)
 +
 +</​code>​
 +
 +
 +<code python | 07/​06-colorize.py>​
 +    """​
 +    Convert an ASCII DEM to an image and colorize
 +    using a heat-map color ramp
 +    """​
 +    import numpy as np
 +    import Image
 +    import ImageOps
 +    import colorsys
 +     
 +    # Source LIDAR DEM file
 +    source = "​lidar.asc"​
 +     
 +    # Output image file
 +    target = "​lidar.bmp"​
 +     
 +    # Load the ASCII DEM into a numpy array
 +    arr = np.loadtxt(source,​ skiprows=6)
 +     
 +    # Convert the numpy array to a PIL image
 +    im = Image.fromarray(arr).convert('​L'​)
 +     
 +    # Enhance the image
 +    im = ImageOps.equalize(im)
 +    im = ImageOps.autocontrast(im)
 +     
 +    # Begin building our color ramp
 +    palette = []
 +     
 +    # Hue, Saturaction,​ Value
 +    # color space
 +    h = .67
 +    s = 1
 +    v = 1
 +     
 +    # We'll step through colors from:
 +    # blue-green-yellow-orange-red.
 +    # Blue=low elevation, Red=high-elevation
 +    step = h/256.0
 +     
 +    # Build the palette
 +    for i in range(256):
 +      rp,gp,bp = colorsys.hsv_to_rgb(h,​s,​v)
 +      r = int(rp*255)
 +      g = int(gp*255)
 +      b = int(bp*255)
 +      palette.extend([r,​g,​b])
 +      h-=step
 +     
 +    # Apply the palette to the image      ​
 +    im.putpalette(palette)
 +     
 +    # Save the image
 +    im.save(target)
 +
 +</​code>​
 +
 +
 +<code python | 07/​07-mesh.py>​
 +    """​
 +    Convert an LAS LIDAR file to a shapefile
 +    by creating a 3D triangle mesh using
 +    Delaunay Triangulation.
 +    """​
 +    # cPickle is used to store
 +    # tessalated triangles
 +    # to save time writing ​
 +    # future shapefiles
 +    import cPickle
 +    import os
 +    import time
 +    import math
 +    # Third-party Python modules:
 +    import numpy as np
 +    import shapefile
 +    from laspy.file import File
 +    import voronoi
 +     
 +    # Source LAS file
 +    source = "​clippedLAS.las"​
 +     
 +    # Output shapefile
 +    target = "​mesh"​
 +     
 +    # Triangles archive
 +    archive = "​triangles.p"​
 +     
 +    # Pyshp archive
 +    pyshp = "​mesh_pyshp.p"​
 +     
 +    # Point class required by
 +    # the voronoi module
 +    class Point:
 +      def __init__(self,​x,​y):​
 +        self.px = x
 +        self.py = y
 +     
 +      def x(self):
 +        return self.px
 +     
 +      def y(self):
 +        return self.py
 +     
 +    # This will be the triangle
 +    # array. ​ Load it from a pickle
 +    # file or use the voronoi module
 +    # to create the triangles.
 +    triangles = None
 +     
 +    if os.path.exists(archive):​
 +      print "​Loading triangle archive..."​
 +      f = open(archive,​ "​rb"​)
 +      triangles = cPickle.load(f)
 +      f.close()
 +      # Open LIDAR LAS file
 +      las = File(source,​ mode="​r"​)  ​
 +    else:
 +      # Open LIDAR LAS file
 +      las = File(source,​ mode="​r"​)  ​
 +      points = []
 +      print "​Assembling points..."  ​
 +      # Pull points from LAS file
 +      for x,y in np.nditer((las.x,​las.y)): ​
 +        points.append(Point(x,​y))
 +      print "​Composing triangles..."​
 +      # Delaunay Triangulation  ​
 +      triangles = voronoi.computeDelaunayTriangulation(points) ​
 +      # Save the triangles to save time if we write more than
 +      # one shapefile.
 +      f = open(archive,​ "​wb"​)
 +      cPickle.dump(triangles,​ f, protocol=2)
 +      f.close()
 +     
 +    print "​Creating shapefile..."​
 +    w = None 
 +    if os.path.exists(pyshp):​
 +      f = open(pyshp, "​rb"​)
 +      w = cPickle.load(f)
 +      f.close()
 +    else:
 +      # PolygonZ shapefile (x,y,z,m)
 +      w = shapefile.Writer(shapefile.POLYGONZ)
 +      w.field("​X1",​ "​C",​ "​40"​)
 +      w.field("​X2",​ "​C",​ "​40"​)
 +      w.field("​X3",​ "​C",​ "​40"​)
 +      w.field("​Y1",​ "​C",​ "​40"​)
 +      w.field("​Y2",​ "​C",​ "​40"​)
 +      w.field("​Y3",​ "​C",​ "​40"​)
 +      w.field("​Z1",​ "​C",​ "​40"​)
 +      w.field("​Z2",​ "​C",​ "​40"​)
 +      w.field("​Z3",​ "​C",​ "​40"​)
 +      tris = len(triangles)
 +      # Loop through shapes and 
 +      # track progress every 10 percent
 +      last_percent = 0
 +      for i in range(tris):​
 +        t = triangles[i]
 +        percent = int((i/​(tris*1.0))*100.0)
 +        if percent % 10.0 == 0 and percent > last_percent:​
 +          last_percent = percent
 +          print "%s %% done - Shape %s/%s at %s" % (percent, i, tris, time.asctime()) ​
 +        part=[]
 +        x1 = las.x[t[0]]
 +        y1 = las.y[t[0]]
 +        z1 = las.z[t[0]]
 +        x2 = las.x[t[1]]
 +        y2 = las.y[t[1]]
 +        z2 = las.z[t[1]]
 +        x3 = las.x[t[2]]
 +        y3 = las.y[t[2]]
 +        z3 = las.z[t[2]]
 +        # Check segments for large triangles
 +        # along the convex hull which is an common
 +        # artificat in Delaunay triangulation
 +        max = 3
 +        if math.sqrt((x2-x1)**2+(y2-y1)**2) > max: continue
 +        if math.sqrt((x3-x2)**2+(y3-y2)**2) > max: continue
 +        if math.sqrt((x3-x1)**2+(y3-y1)**2) > max: continue
 +        part.append([x1,​y1,​z1,​0])
 +        part.append([x2,​y2,​z2,​0])
 +        part.append([x3,​y3,​z3,​0])
 +        w.poly(parts=[part])
 +        w.record(x1,​x2,​x3,​y1,​y2,​y3,​z1,​z2,​z3)
 +      print "​Saving shapefile..."​
 +      # Pickle the Writer in case something
 +      # goes wrong. Be sure to delete this
 +      # file to recreate teh shapefile.
 +      f = open(pyshp, "​wb"​)
 +      cPickle.dump(w,​ f, protocol=2)
 +      f.close()
 +    w.save(target)
 +    print "​Done."​
 +
 +
 +</​code>​
 +
 +===== More advanced modules =====
 +
 +
 +<code python | 08/​01-ndvi.py>​
 +    """​
 +    Output a normalized vegetative index
 +    """​
 +    import gdal, gdalnumeric,​ ogr
 +    import Image, ImageDraw
 +     
 +    def imageToArray(i):​
 +      """​
 +      Converts a Python Imaging Library ​
 +      array to a gdalnumeric image.
 +      """​
 +      a=gdalnumeric.numpy.fromstring(i.tostring(),'​b'​)
 +      a.shape=i.im.size[1],​ i.im.size[0]
 +      return a
 +     
 +    def world2Pixel(geoMatrix,​ x, y):
 +      """​
 +      Uses a gdal geomatrix (gdal.GetGeoTransform()) ​
 +      to calculate the pixel location of a 
 +      geospatial coordinate ​
 +      """​
 +      ulX = geoMatrix[0]
 +      ulY = geoMatrix[3]
 +      xDist = geoMatrix[1]
 +      yDist = geoMatrix[5]
 +      rtnX = geoMatrix[2]
 +      rtnY = geoMatrix[4]
 +      pixel = int((x - ulX) / xDist)
 +      line = int((ulY - y) / xDist)
 +      return (pixel, line)    ​
 +     
 +    # Multispectral image used 
 +    # to create the NDVI. Must
 +    # have red and infrared
 +    # bands
 +    source = "​farm.tif"​
 +     
 +    # Output geotiff file name
 +    target = "​ndvi.tif"​
 +     
 +    # Load the source data as a gdalnumeric array
 +    srcArray = gdalnumeric.LoadFile(source)
 +     
 +    # Also load as a gdal image to 
 +    # get geotransform (world file) info
 +    srcImage = gdal.Open(source)
 +    geoTrans = srcImage.GetGeoTransform()
 +     
 +    # Red and infrared (or near infrared) bands
 +    r = srcArray[1]
 +    ir = srcArray[2]
 +     
 +    ## Clip a field out of the bands using a
 +    ## field boundary shapefile
 +     
 +    # Create an OGR layer from a Field boundary shapefile
 +    field = ogr.Open("​field.shp"​)
 +    # Must define a "​layer"​ to keep OGR happy
 +    lyr = field.GetLayer("​field"​)
 +    # Only one polygon in this shapefile
 +    poly = lyr.GetNextFeature()
 +     
 +    # Convert the layer extent to image pixel coordinates
 +    minX, maxX, minY, maxY = lyr.GetExtent()
 +    ulX, ulY = world2Pixel(geoTrans,​ minX, maxY)
 +    lrX, lrY = world2Pixel(geoTrans,​ maxX, minY)
 +     
 +    # Calculate the pixel size of the new image
 +    pxWidth = int(lrX - ulX)
 +    pxHeight = int(lrY - ulY)
 +     
 +    # Create a blank image of the correct size
 +    # that will serve as our mask
 +    clipped = gdalnumeric.numpy.zeros((3,​ pxHeight, pxWidth), \
 +    gdalnumeric.numpy.uint8)
 +    #mmask = gdalnumeric.zeros((3,​ pxHeight, pxWidth), gdalnumeric.UnsignedInt8)
 +    #rgb = rgb.astype(gdalnumeric.UnsignedInt8)
 +    rClip = r[ulY:lrY, ulX:lrX]
 +    irClip = ir[ulY:lrY, ulX:lrX]
 +     
 +    # Create a new geomatrix for the image
 +    geoTrans = list(geoTrans)
 +    geoTrans[0] = minX
 +    geoTrans[3] = maxY
 +     
 +    # Map points to pixels for drawing ​
 +    # the field boundary on a blank
 +    # 8-bit, black and white, mask image.
 +    points = []
 +    pixels = []
 +    # Grab the polygon geometry
 +    geom = poly.GetGeometryRef()
 +    pts = geom.GetGeometryRef(0)
 +    # Loop through geometry and turn
 +    # the points into an easy-to-manage ​
 +    # Python list 
 +    for p in range(pts.GetPointCount()):​
 +      points.append((pts.GetX(p),​ pts.GetY(p)))
 +    # Loop through the points and map to pixels.
 +    # Append the pixels to a pixel list
 +    for p in points:
 +      pixels.append(world2Pixel(geoTrans,​ p[0], p[1]))
 +    # Create the raster polygon image
 +    rasterPoly = Image.new("​L",​ (pxWidth, pxHeight), 1)
 +    # Create a PIL drawing object
 +    rasterize = ImageDraw.Draw(rasterPoly)
 +    # Dump the pixels to the image
 +    rasterize.polygon(pixels,​ 0)
 +    # Hand the image back to gdal/​gdalnumeric
 +    # so we can use it as an array mask
 +    mask = imageToArray(rasterPoly) ​   ​
 +    # Clip the red band using the mask   
 +    rClip = gdalnumeric.numpy.choose(mask,​ \
 +      (rClip, 0)).astype(gdalnumeric.numpy.uint8)
 +    # Clip the infrared band using the mask
 +    irClip = gdalnumeric.numpy.choose(mask,​ \
 +      (irClip, 0)).astype(gdalnumeric.numpy.uint8)
 +     
 +    # We don't care about numpy warnings
 +    # due to NaN values from clipping
 +    gdalnumeric.numpy.seterr(all="​ignore"​)
 +     
 +    # NDVI equation: (infrared - red) / (infrared + red)
 +    # *1.0 converts values to floats, ​
 +    # +1.0 prevents ZeroDivisionErrors ​
 +    ndvi = 1.0 * (irClip - rClip) / irClip + rClip + 1.0
 +     
 +    # Remove any NaN values from the final product
 +    ndvi = gdalnumeric.numpy.nan_to_num(ndvi)
 +     
 +    # Save ndvi as tiff
 +    gdalnumeric.SaveArray(ndvi,​ target, \
 +      format="​GTiff",​ prototype=source)
 +
 +
 +</​code>​
 +
 +
 +<code python | 08/​02-ndvi-classify.py>​
 +    """​
 +    Classify an NDVI tiff using 7 classes ​
 +    by "​pushing"​ the NDVI through
 +    masks defined by the desired ​
 +    range of values for each class.
 +    """​
 +    import gdalnumeric as gd
 +    import operator
 +     
 +    def histogram(a,​ bins=range(0,​256)):​
 +      """​
 +      Histogram function for multi-dimensional array.
 +      a = array
 +      bins = range of numbers to match 
 +      """​
 +      fa = a.flat
 +      n = gd.numpy.searchsorted(gd.numpy.sort(fa),​ bins)
 +      n = gd.numpy.concatenate([n,​ [len(fa)]])
 +      hist = n[1:​]-n[:​-1] ​
 +      return hist
 +     
 +    def stretch(a):
 +      """​
 +      Performs a histogram stretch on a gdalnumeric array image.
 +      """​
 +      hist = histogram(a)
 +      lut = []
 +      for b in range(0, len(hist), 256):
 +        # step size
 +        step = reduce(operator.add,​ hist[b:​b+256]) / 255
 +        # create equalization lookup table
 +        n = 0
 +        for i in range(256):
 +          lut.append(n / step)
 +          n = n + hist[i+b]
 +      gd.numpy.take(lut,​ a, out=a)
 +      return a
 +     
 +    # NDVI output from ndvi script
 +    source = "​ndvi.tif"​
 +    # Target file name for classified
 +    # image image
 +    target = "​ndvi_color.tif"​
 +     
 +    # Load the image into an array
 +    ndvi = gd.LoadFile(source).astype(gd.numpy.uint8)
 +    # Peform a histogram stretch so we are able to
 +    # use all of the classes
 +    ndvi = stretch(ndvi)
 +     
 +    # Create a blank 3-band image the same size as the ndvi
 +    rgb = gd.numpy.zeros((3,​ len(ndvi), len(ndvi[0])),​ gd.numpy.uint8)
 +     
 +    # Class list with ndvi upper range values.
 +    # Note the lower and upper values are listed on the ends
 +    classes = [58,​73,​110,​147,​184,​220,​255]
 +     
 +    # Color look-up table (lut)
 +    # The lut must match the number of classes
 +    # Specified as R,G,B tuples from dark brown to dark green
 +    lut = [[120,​69,​25],​ [255,​178,​74],​ [255,​237,​166],​ [173,​232,​94],​
 +           ​[135,​181,​64],​ [3,156,0], [1,100,0]]
 +     
 +    # Starting value of the first class
 +    start = 1
 +     
 +    # Process all classes.
 +    for i in range(len(classes)):​
 +        mask = gd.numpy.logical_and(\
 +        start <= ndvi, ndvi <= classes[i])
 +        for j in range(len(lut[i])):​
 +            rgb[j] = gd.numpy.choose(mask,​ \
 +              (rgb[j], lut[i][j]))
 +        start = classes[i]+1 ​   ​
 +     
 +    # Save a geotiff image of the colorized ndvi.
 +    gd.SaveArray(rgb,​ target, format="​GTiff",​ prototype=source)
 +
 +
 +</​code>​
 +
 +
 +<code python | 08/​03-flood-fill.py>​
 +    """​
 +    Crawls a terrain raster from a starting
 +    point and "​floods"​ everything at the same
 +    or lower elevation by producing a mask
 +    image of 1,0 values.
 +    """​
 +     
 +    import numpy as np
 +    from linecache import getline
 +     
 +    def floodFill(c,​r,​mask):​
 +      """​
 +      Crawls a mask array containing
 +      only 1 and 0 values from the
 +      starting point (c=column,
 +      r=row - a.k.a. x,y) and returns
 +      an array with all 1 values
 +      connected to the starting cell.
 +      This algorithm performs a 4-way
 +      check non-recursively.  ​
 +      """​
 +      # cells already filled
 +      filled = set()
 +      # cells to fill
 +      fill = set()
 +      fill.add((c,​r))
 +      width = mask.shape[1]-1
 +      height = mask.shape[0]-1
 +      # Our output inundation array
 +      flood = np.zeros_like(mask,​ dtype=np.int8)
 +      # Loop through and modify the cells which
 +      # need to be checked.
 +      while fill:
 +        # Grab a cell
 +        x,y = fill.pop()
 +        if y == height or x == width or x < 0 or y < 0:
 +          # Don't fill
 +          continue
 +        if mask[y][x] == 1:
 +          # Do fill
 +          flood[y][x]=1
 +          filled.add((x,​y))
 +          # Check neighbors for 1 values
 +          west =(x-1,y)
 +          east = (x+1,y)
 +          north = (x,y-1)
 +          south = (x,y+1)
 +          if not west in filled:
 +            fill.add(west)
 +          if not east in filled: ​     ​
 +            fill.add(east)
 +          if not north in filled: ​
 +            fill.add(north)
 +          if not south in filled: ​
 +            fill.add(south)
 +      return flood 
 +     
 +    source = "​terrain.asc"​
 +    target = "​flood.asc"​
 +     
 +    print "​Opening image..."​
 +    img = np.loadtxt(source,​ skiprows=6)
 +    print "Image opened"​
 +     
 +    a = np.where(img<​70,​ 1,0)
 +    print "Image masked"​
 +     
 +    # Parse the headr using a loop and
 +    # the built-in linecache module
 +    hdr = [getline(source,​ i) for i in range(1,7)]
 +    values = [float(h.split("​ "​)[-1].strip()) for h in hdr]
 +    cols,​rows,​lx,​ly,​cell,​nd = values
 +    xres = cell
 +    yres = cell * -1
 +     
 +    # Starting point for the 
 +    # flood inundation ​
 +    sx = 2582
 +    sy = 2057
 +     
 +    print "​Beginning flood fill" ​     ​
 +    fld = floodFill(sx,​sy,​ a)
 +    print "​Finished Flood fill"
 +     
 +    header=""​
 +    for i in range(6):
 +      header += hdr[i]
 +     
 +    print "​Saving grid"
 +    # Open the output file, add the hdr, save the array
 +    with open(target,​ "​wb"​) as f:
 +      f.write(header)
 +      np.savetxt(f,​ fld, fmt="​%1i"​)
 +    print "​Done!"​
 +
 +
 +</​code>​
 +
 +
 +<code python | 08/​04-least-cost_model.py>​
 +"""​
 +Provides a simple command-line-output
 +version of the least cost path solution
 +using randomly-generated notional arrays.
 +"""​
 + 
 +import numpy as np
 + 
 +# Our A* search algorithm  ​
 +def astar(start,​ end, h, g):
 +  cset = set()
 +  oset = set()
 +  path = set()
 +  oset.add(start)
 +  while oset:
 +    cur = oset.pop()
 +    if cur == end:
 +      return path
 +    cset.add(cur)
 +    path.add(cur)
 +    options = []
 +    y1 = cur[0]
 +    x1 = cur[1]
 +    if y1 > 0: 
 +      options.append((y1-1,​ x1))
 +    if y1 < h.shape[0]-1: ​
 +      options.append((y1+1,​ x1))
 +    if x1 > 0: 
 +      options.append((y1,​ x1-1))
 +    if x1 < h.shape[1]-1: ​
 +      options.append((y1,​ x1+1))
 +    if end in options:
 +      return path
 +    best = options[0]
 +    cset.add(options[0])
 +    for i in range(1,​len(options)):​
 +      option = options[i]
 +      if option in cset:
 +        continue
 +      elif h[option] <= h[best]:
 +        best = option ​  
 +        cset.add(option)
 +      elif g[option] < g[best]:
 +        best = option
 +        cset.add(option)
 +      else:
 +        cset.add(option)
 +    print best, ", ", h[best], ", ", g[best]
 +    oset.add(best)
 +  return []
 + 
 + 
 +# Width and height
 +# of grids
 +w = 5
 +h = 5
 + 
 +# Start location:
 +# Lower left of grid
 +start = (h-1, 0)
 + 
 +# End location:
 +# Top right of grid
 +dx = w-1
 +dy = 0
 + 
 +# Blank grid
 +a = np.zeros((w,​h))
 + 
 +# Distance grid
 +dist = np.zeros(a.shape,​ dtype=np.int8)
 + 
 +# Calculate distance for all cells
 +for y,x in np.ndindex(a.shape):​
 +  dist[y][x] = abs((dx-x)+(dy-y))
 + 
 +# "​Terrain"​ is a random value between 1-16.
 +# Add to the distance grid to calculate
 +# The cost of moving to a cell
 +cost = np.random.randint(1,​16,​(w,​h)) + dist
 + 
 +print "COST GRID (Value + Distance)"​
 +print cost
 +print
 + 
 +print "​(Y,​X),​ HEURISTIC, DISTANCE"​
 +# Find the path
 +path = astar(start,​(dy,​dx),​cost,​ dist)
 +print
 + 
 +# Create and populate the path grid
 +path_grid = np.zeros(cost.shape,​ dtype=np.uint8)
 +for y,x in path:
 +  path_grid[y][x]=1
 +path_grid[dy][dx]=1
 + 
 +print "PATH GRID: 1=path"​
 +print path_grid
 +</​code>​
 +
 +
 +<code python | 08/​05-least-cost-path.py>​
 +    """​
 +    Calculates the least cost path
 +    over a terrain grid and outputs
 +    another raster of 1,0 values
 +    defining the path.
 +    """​
 +     
 +    import numpy as np
 +    import math
 +    from linecache import getline
 +     
 +    def e_dist(p1,​p2):​
 +      """​
 +      Takes two points and returns
 +      the euclidian distance
 +      """​
 +      x1,y1=p1
 +      x2,y2=p2
 +      distance = math.sqrt((x1-x2)**2+(y1-y2)**2)
 +      return int(distance)
 +     
 +    def weighted_score(cur,​ node, h, start, end):
 +      """​
 +      Provides a weighted score by comparing the
 +      current node with a neighboring node. Loosely
 +      based on on the Nisson score concept: f=g+h
 +      In this case, the "​h"​ value, or "​heuristic",​
 +      is the elevation value of each node.     
 +      """​
 +      score = 0
 +      # current node elevation
 +      cur_h = h[cur]
 +      # current node distance from end
 +      cur_g = e_dist(cur,​end)
 +      # current node distance from 
 +      cur_d = e_dist(cur,​start)
 +      # neighbor node elevation
 +      node_h = h[node]
 +      # neighbor node distance from end
 +      node_g = e_dist(node,​end)
 +      # neighbor node distance from start
 +      node_d = e_dist(node,​ start)
 +      # Compare values with the heighest
 +      # weight given to terrain followed
 +      # by progress towards the goal.
 +      if node_h < cur_h:
 +        score += cur_h-node_h
 +      if node_g < cur_g:
 +        score += 10
 +      if node_d > cur_d:
 +        score += 10
 +      return score
 +     
 +    def astar(start,​ end, h):
 +      """​
 +      A-Star (or A*) search algorithm.
 +      Moves through nodes in a network
 +      (or grid), scores each node'​s ​
 +      neighbors, and goes to the node
 +      with the best score until it finds
 +      the end.  A* is an evolved Dijkstra
 +      algorithm.
 +      """​
 +      # Closed set of nodes to avoid
 +      cset = set()
 +      # Open set of nodes to evaluate
 +      oset = set()
 +      # Output set of path nodes
 +      path = set()
 +      # Add the starting point to
 +      # to begin processing
 +      oset.add(start)
 +      while oset:
 +        # Grab the next node
 +        cur = oset.pop()
 +        # Return if we're at the end
 +        if cur == end:
 +          return path
 +        # Close off this node to future
 +        # processing
 +        cset.add(cur)
 +        # The current node is always
 +        # a path node by definition
 +        path.add(cur)
 +        # List to hold neighboring
 +        # nodes for processing
 +        options = []
 +        # Grab all of the neighbors
 +        y1 = cur[0]
 +        x1 = cur[1]
 +        if y1 > 0: 
 +          options.append((y1-1,​ x1))
 +        if y1 < h.shape[0]-1: ​
 +          options.append((y1+1,​ x1))
 +        if x1 > 0: 
 +          options.append((y1,​ x1-1))
 +        if x1 < h.shape[1]-1: ​
 +          options.append((y1,​ x1+1))
 +        if x1 > 0 and y1 > 0:
 +          options.append((y1-1,​ x1-1))
 +        if y1 < h.shape[0]-1 and x1 <  h.shape[1]-1:​
 +          options.append((y1+1,​ x1+1))
 +        if y1 < h.shape[0]-1 and x1 > 0:
 +          options.append((y1+1,​ x1-1))
 +        if y1 > 0 and x1 <  h.shape[1]-1:​
 +          options.append((y1-1,​ x1+1))  ​
 +        # If the end is a neighbor, return
 +        if end in options:
 +          return path
 +        # Store the best known node
 +        best = options[0]
 +        # Begin scoring neighbors
 +        best_score = weighted_score(cur,​ best, h, start, end)    ​
 +        # process the other 7 neighbors
 +        for i in range(1,​len(options)):​
 +          option = options[i]
 +          # Make sure the node is new
 +          if option in cset:
 +            continue
 +          else:
 +            # Score the option and compare it to the best known
 +            option_score = weighted_score(cur,​ option, h, start, end)
 +            if option_score > best_score:
 +              best = option
 +              best_score = option_score
 +            else:
 +              # If the node isn't better seal it off
 +              cset.add(option)
 +        # Uncomment this print statement to watch
 +        # the path develop in real time:
 +        # print best, e_dist(best,​end)
 +     
 +        # Add the best node to the open set
 +        oset.add(best)
 +      return []
 +     
 +    # Our terrain data
 +    source = "​dem.asc"​
 +     
 +    # Output file name
 +    # for the path raster
 +    target = "​path.asc"​
 +     
 +    print "​Opening %s..." % source
 +    cost = np.loadtxt(source,​ skiprows=6)
 +    print "​Opened %s." % source
 +     
 +    # Parse the header
 +    hdr = [getline(source,​ i) for i in range(1,7)]
 +    values = [float(ln.split("​ "​)[-1].strip()) for ln in hdr]
 +    cols,​rows,​lx,​ly,​cell,​nd = values
 +     
 +    # Starting column, row
 +    sx = 1006
 +    sy = 954
 +     
 +    # Ending column, row
 +    dx = 303 
 +    dy = 109
 +     
 +    print "​Searching for path..."​
 +    p = astar((sy,​sx),​(dy,​dx),​cost)
 +    print "Path found."​
 +    print "​Creating path grid..."​
 +    path = np.zeros(cost.shape)
 +    print "​Plotting path..."​
 +    for y,x in p:
 +      path[y][x]=1
 +    path[dy][dx]=1
 +    print "Path plotted."​
 +     
 +    print "​Saving %s..." % target ​   ​
 +    header=""​
 +    for i in range(6):
 +      header += hdr[i]
 +     
 +    # Open the output file, add the hdr, save the array
 +    with open(target,​ "​wb"​) as f:
 +      f.write(header)
 +      np.savetxt(f,​ path, fmt="​%4i"​)
 +     
 +    print "​Done!"​
 +
 +
 +</​code>​
 +
 +<code python | 09/​01-nextbus.py>​
 +import urllib
 +from xml.dom import minidom
 +
 +# Nextbus API command mode
 +command = "​vehicleLocations"​
 +
 +# Nextbus customer to query
 +agency = "​chapel-hill"​
 +
 +# Bus we want to query
 +route = "​A"​
 +
 +# Time in milliseconds since the 
 +# 1970 epoch time.  All tracks
 +# after this time will be returned.
 +# 0 only returns data for the last
 +# 15 minutes
 +epoch = "​0"​
 +
 +## Build our query url
 +#
 +# webservices base url
 +url = "​http://​webservices.nextbus.com"​
 +# web service path
 +url += "/​service/​publicXMLFeed?"​
 +# service command/​mode
 +url += "​command="​ + command
 +# agency ​
 +url += "&​a="​ + agency
 +url +=  "&​r="​ + route
 +url += "&​t="​ + epoch
 +
 +# Access the REST URL
 +feed = urllib.urlopen(url)
 +
 +
 +if feed:
 +  # Parse the xml feed
 +  xml = minidom.parse(feed)
 +  # Get the vehicle tags
 +  vehicles = xml.getElementsByTagName("​vehicle"​)
 +  # Get the most recent one. Normally there will
 +  # be only one.
 +  if vehicles:
 +    bus = vehicles.pop()
 +    # Print the bus latitude and longitude
 +    att = bus.attributes
 +    print att["​lon"​].value,​ ",",​ att["​lat"​].value
 +  else:
 +    print "No vehicles found." ​
 +
 +</​code>​
 +
 +<code python | 09/​02-nextma>​
 +import urllib
 +from xml.dom import minidom
 +import time
 +
 +def nextbus(a, r, c="​vehicleLocations",​ e=0):
 +  """​Returns the most recent latitude and 
 +  longitude of the selected bus line using
 +  the NextBus API (nbapi)"""​
 +  nbapi = "​http://​webservices.nextbus.com"​
 +  nbapi += "/​service/​publicXMLFeed?"​
 +  nbapi += "​command=%s&​a=%s&​r=%s&​t=%s"​ % (c,a,r,e)
 +  xml = minidom.parse(urllib.urlopen(nbapi))
 +  bus=xml.getElementsByTagName("​vehicle"​)
 +  if bus:    ​
 +    at = bus[0].attributes
 +    return(at["​lat"​].value,​ at["​lon"​].value)
 +  else: return (False, False)
 +
 +def nextmap(a, r, mapimg, index):
 +  """​Plots a nextbus location on a map image
 +  and saves it to disk using the OpenStreetMap
 +  Static Map API (osmapi)"""​
 +  # Fetch the latest bus location
 +  lat, lon = nextbus(a, r)
 +  if not lat:
 +    return False
 +  # Base url + service path
 +  osmapi = "​http://​staticmap.openstreetmap.de/"​
 +  osmapi += "​staticmap.php?​maptype=mapnik&"​
 +  # Map Image width and height in pixels
 +  osmapi += "​size=800x600"​ + "&"​
 +  # Center the map on the bus location
 +  osmapi += "​center=%s,​%s"​ % (lat,lon) + "&"​
 +  # Map zoom level (between 1-18)
 +  osmapi += "​zoom=16"​ + "&"​
 +  # Bus mark location ​
 +  osmapi += "​markers=%s,​%s,​lightblue%d"​ % (lat, lon, index) ​
 +  print osmapi
 +  img = urllib.urlopen(osmapi)
 +  # Save the map image
 +  with open(mapimg + "​%d.png"​ % index, "​wb"​) as f:
 +    f.write(img.read())
 +  return True
 +
 +# Nextbus API agency and bus line variables
 +agency = "​chapel-hill"​
 +route = "​A"​
 +
 +# Name of map image to save as PNG
 +nextimg = "​nextmap"​
 +
 +# Number of updates we want to make
 +requests = 3
 +
 +# How often we want to update (seconds)
 +freq = 5
 +
 +# Map the bus location every few seconds ​
 +for i in range(requests):​
 +  success = nextmap(agency,​ route, nextimg, i+1)
 +  if not success:
 +    print "No data available."​
 +    continue
 +  print "Saved map %s at %s" % (i, time.asctime())
 +  time.sleep(freq)
 +</​code>​
  
wiki/python/basicpython3.txt ยท Last modified: 2018/06/01 20:13 (external edit)