Python：利用matplotlib-natgrid工具箱对大量点进行快速插值

Question

我正在使用matplotlib natgrid工具箱来互操作x、y、z点。我的数据集超过500万点。我用一小块区域（大约900.000,00点）尝试了我的代码。使用natgrid的时间是44分钟。在

有人知道一种提高速度的方法，还是另一种在时间上更有效的方法？对于二维插值，有太多的数据点需要插值

提前感谢您的帮助和建议

import shapefile #
import os #
import glob #
import math #
import numpy #
import numpy as np #
import matplotlib.nxutils as nx #
import collections
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
import matplotlib.delaunay
from liblas import file as lasfile #
from shapely.geometry import Polygon #
from osgeo import gdal, osr, ogr #
from gdalconst import * #
from matplotlib.mlab import griddata
from collections import OrderedDict


    def LAS2DTM(inFile,outFile,gridSize=1,dtype="GDT_Float32",nodata=-9999.00,BBOX=None,EPSG=None):
        if BBOX == None:
            X = []
            Y = []
            for p in lasfile.File(inFile,None,'r'):
                X.append(p.x)
                Y.append(p.y)
            xmax, xmin = max(X),min(X)
            ymax, ymin = max(Y), min(Y)
            del X,Y
        else:
            xmax,xmin,ymax,ymin = BBOX[0],BBOX[1],BBOX[2],BBOX[3]
        # number of row and columns
        nx = int(math.ceil(abs(xmax - xmin)/gridSize))
        ny = int(math.ceil(abs(ymax - ymin)/gridSize))
        # Create an array to hold the number of points in each pixel
        cnts = np.zeros((ny, nx))
        # Create an array to hold the values
        data = np.zeros((ny, nx))
        # read all points
        x = []
        y = []
        z = []
        for p in lasfile.File(inFile,None,'r'):
            x.append(p.x)
            y.append(p.y)
            z.append(p.z)
            # Compute the x and y offsets for where this point would be in the raster
            dx = int((p.x - xmin)/gridSize)
            dy = int((ymax - p.y)/gridSize)
            # Add the z value to the total for that pixel
            data[dy,dx] += p.z
            # Add 1 to our count for that pixel
            cnts[dy,dx] += 1
        # ingore Error message
        np.seterr(invalid='ignore')
        # Compute the averages
        data = data/cnts
        del cnts
        # remove all duplicate points from a X,Y,Z file that have identical x and y coordinates
        # The first point survives, all subsequent duplicates are removed.
        tmp = OrderedDict()
        for point in zip(x, y, z):
           a = tmp.setdefault(point[:2], point)
        mypoints = tmp.values()
        del x,y,z
        points_zipped = zip(*mypoints)
        del mypoints
        xvals = np.array(points_zipped[0])
        yvals = np.array(points_zipped[1])
        zvals = np.array(points_zipped[2])
        del points_zipped
        # define grid.
        xi = np.linspace(xmin, xmax, nx)
        yi = np.linspace(ymin, ymax, ny)
        # create a meshgrid
        xi, yi = np.meshgrid(xi, yi)
        # grid the data.
        zi = griddata(xvals,yvals,zvals,xi,yi,interp='nn')
        # convert "numpy.ma.core.MaskedArray" in a "np.array"
        zi = np.array(zi)
        # mask a numpy.ndarray with another numpy.ndarray
        data[np.isnan(data)] = zi[np.isnan(data)]
        # Create gtif
        if dtype == "GDT_Unknown": # Unknown or unspecified type
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_Unknown)
        elif dtype == "GDT_Byte": # Eight bit unsigned integer
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_Byte)
        elif dtype == "GDT_UInt16": # Sixteen bit unsigned integer
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_UInt16)
        elif dtype == "GDT_Int16": # Sixteen bit signed integer
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_Int16)
        elif dtype == "GDT_UInt32": # Thirty two bit unsigned integer
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_UInt32)
        elif dtype == "GDT_Int32": # Thirty two bit signed integer
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_Int32)
        elif dtype == "GDT_Float32": # Thirty two bit floating point
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_Float32)
        elif dtype == "GDT_Float64": # Sixty four bit floating point
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_Float64)
        elif dtype == "GDT_CInt16": # Complex Int16
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_CInt16)
        elif dtype == "GDT_CInt32": # Complex Int32
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_CInt32)
        elif dtype == "GDT_CFloat32": # Complex Float32
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_CFloat32)
        elif dtype == "GDT_CFloat64": # Complex Float64
            target_ds = gdal.GetDriverByName('GTiff').Create(outFile, nx,ny, 1, gdal.GDT_CFloat64)
        # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
        target_ds.SetGeoTransform((xmin, gridSize, 0,ymax, 0, -gridSize))
        # set the reference info
        if EPSG is None:
            # Source has no projection (needs GDAL >= 1.7.0 to work)
            target_ds.SetProjection('LOCAL_CS["arbitrary"]')
        else:
            proj = osr.SpatialReference()
            proj.ImportFromEPSG(EPSG)
            # Make the target raster have the same projection as the source
            target_ds.SetProjection(proj.ExportToWkt())
        # write the band
        target_ds.GetRasterBand(1).WriteArray(data)
        target_ds.GetRasterBand(1).SetNoDataValue(nodata)
        target_ds = None

Answer 1

尝试不同类型的二维插值，例如可以尝试interp2d。在