Python osgeo.ogr 模块，CreateGeometryFromWkt() 实例源码

def get_modis_tile_list(ds):
    """Helper function to identify MODIS tiles that intersect input geometry

    modis_gird.py contains dictionary of tile boundaries (tile name and WKT polygon ring from bbox)

    See: https://modis-land.gsfc.nasa.gov/MODLAND_grid.html
    """
    from demcoreg import modis_grid
    modis_dict = modis_grid.modis_dict
    for key in modis_dict:
        modis_dict[key] = ogr.CreateGeometryFromWkt(modis_dict[key])
    geom = geolib.ds_geom(ds)
    geom_dup = geolib.geom_dup(geom)
    ct = osr.CoordinateTransformation(geom_dup.GetSpatialReference(), geolib.wgs_srs)
    geom_dup.Transform(ct)
    tile_list = []
    for key, val in list(modis_dict.items()):
        if geom_dup.Intersects(val):
            tile_list.append(key)
    return tile_list

def shp2geom(shp_fn):
    """Extract geometries from input shapefile

    Need to handle multi-part geom: http://osgeo-org.1560.x6.nabble.com/Multipart-to-singlepart-td3746767.html
    """
    ds = ogr.Open(shp_fn)
    lyr = ds.GetLayer()
    srs = lyr.GetSpatialRef()
    lyr.ResetReading()
    geom_list = []
    for feat in lyr:
        geom = feat.GetGeometryRef()
        geom.AssignSpatialReference(srs)
        #Duplicate the geometry, or segfault
        #See: http://trac.osgeo.org/gdal/wiki/PythonGotchas
        #g = ogr.CreateGeometryFromWkt(geom.ExportToWkt())
        #g.AssignSpatialReference(srs)
        g = geom_dup(geom)
        geom_list.append(g)
    #geom = ogr.ForceToPolygon(' '.join(geom_list))    
    #Dissolve should convert multipolygon to single polygon 
    #return geom_list[0]
    ds = None
    return geom_list

def ds_geom(ds, t_srs=None):
    """Return dataset bbox envelope as geom
    """
    gt = ds.GetGeoTransform()
    ds_srs = get_ds_srs(ds)
    if t_srs is None:
        t_srs = ds_srs
    ns = ds.RasterXSize
    nl = ds.RasterYSize
    x = np.array([0, ns, ns, 0, 0], dtype=float)
    y = np.array([0, 0, nl, nl, 0], dtype=float)
    #Note: pixelToMap adds 0.5 to input coords, need to account for this here
    x -= 0.5
    y -= 0.5
    mx, my = pixelToMap(x, y, gt)
    geom_wkt = 'POLYGON(({0}))'.format(', '.join(['{0} {1}'.format(*a) for a in zip(mx,my)]))
    geom = ogr.CreateGeometryFromWkt(geom_wkt)
    geom.AssignSpatialReference(ds_srs)
    if not ds_srs.IsSame(t_srs):
        geom_transform(geom, t_srs)
    return geom

def get_features_as_geojson(layer, bbox=None, union=False):
    """ Get features in this layer and return as GeoJSON """
    if bbox is not None:
        layer.SetSpatialFilterRect(bbox[0], bbox[3], bbox[2], bbox[1])
    poly = ogr.Geometry(ogr.wkbPolygon)
    if union:
        for feature in layer:
            geom = feature.GetGeometryRef()
            #  required for ogr2
            if hasattr(geom, 'GetLinearGeometry'):
                geom = geom.GetLinearGeometry()
            poly = poly.Union(geom)
    if bbox is not None:
        wkt = "POLYGON ((%s %s, %s %s, %s %s, %s %s, %s %s))" % \
              (bbox[0], bbox[1], bbox[2], bbox[1], bbox[2], bbox[3], bbox[0], bbox[3], bbox[0], bbox[1])
        bbox_wkt = ogr.CreateGeometryFromWkt(wkt)
        poly = poly.Intersection(bbox_wkt)
    return json.loads(poly.ExportToJson())

def get_idx_as_shp(self, path, gt, wkt):
        '''
        Exports a Shapefile containing the locations of the extracted
        endmembers. Assumes the coordinates are in decimal degrees.
        '''
        coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)

        driver = ogr.GetDriverByName('ESRI Shapefile')
        ds = driver.CreateDataSource(path)
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)

        layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint)
        for pair in coords:
            feature = ogr.Feature(layer.GetLayerDefn())

            # Create the point from the Well Known Text
            point = ogr.CreateGeometryFromWkt('POINT(%f %f)' %  pair)
            feature.SetGeometry(point) # Set the feature geometry
            layer.CreateFeature(feature) # Create the feature in the layer
            feature.Destroy() # Destroy the feature to free resources

        # Destroy the data source to free resources
        ds.Destroy()

def readwktcsv(csv_path,removeNoBuildings=True, groundTruthFile=True):
    #
    # csv Format Expected = ['ImageId', 'BuildingId', 'PolygonWKT_Pix', 'PolygonWKT_Geo']
    # returns list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'poly': poly}
    # image_id is a string,
    # BuildingId is an integer,
    # poly is a ogr.Geometry Polygon

    buildinglist = []
    with open(csv_path, 'rb') as csvfile:
        building_reader = csv.reader(csvfile, delimiter=',', quotechar='"')
        next(building_reader, None)  # skip the headers
        for row in building_reader:

            if removeNoBuildings:
                if int(row[1]) != -1:
                    polyPix = ogr.CreateGeometryFromWkt(row[2])
                    if groundTruthFile:
                        polyGeo = ogr.CreateGeometryFromWkt(row[3])
                    else:
                        polyGeo = []
                    buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
                                         'polyGeo': polyGeo,
                                         })

            else:

                polyPix = ogr.CreateGeometryFromWkt(row[2])
                if groundTruthFile:
                    polyGeo = ogr.CreateGeometryFromWkt(row[3])
                else:
                    polyGeo = []
                buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
                                     'polyGeo': polyGeo,
                                     })

    return buildinglist

def pixelWKTToGeoWKT(geomWKT, inputRaster, targetSR='', geomTransform='', breakMultiPolygonPix=False):
    # Returns  GeoCoordinateList from PixelCoordinateList



    geomPix = ogr.CreateGeometryFromWkt(geomWKT)
    geomGeoList = pixelGeomToGeoGeom(geomPix, inputRaster, targetSR=targetSR,
                                 geomTransform=geomTransform, breakMultiPolygonPix=breakMultiPolygonPix)

    return geomGeoList

def geom_dup(geom):
    """Create duplicate geometry

    Needed to avoid segfault when passing geom around. See: http://trac.osgeo.org/gdal/wiki/PythonGotchas
    """
    g = ogr.CreateGeometryFromWkt(geom.ExportToWkt())
    g.AssignSpatialReference(geom.GetSpatialReference())
    return g 

#This should be a function of a new geom class
#Assumes geom has srs defined
#Modifies geom in place

def bbox2geom(bbox, t_srs=None):
    #Check bbox
    #bbox = numpy.array([-180, 180, 60, 90])
    x = [bbox[0], bbox[0], bbox[1], bbox[1], bbox[0]]
    y = [bbox[2], bbox[3], bbox[3], bbox[2], bbox[2]]
    geom_wkt = 'POLYGON(({0}))'.format(', '.join(['{0} {1}'.format(*a) for a in zip(x,y)]))
    geom = ogr.CreateGeometryFromWkt(geom_wkt)
    if t_srs is not None and not wgs_srs.IsSame(t_srs):
        ct = osr.CoordinateTransformation(t_srs, wgs_srs)
        geom.Transform(ct)
        geom.AssignSpatialReference(t_srs)
    return geom

def intersection(self,filename):
        """ Return coordinates of intersection between this image and another.

        :param filename : path to the second image (or another GeoRaster instance)
        :type filename: str, georaster.__Raster

        :returns: extent of the intersection between the 2 images \
        (xmin, xmax, ymin, ymax)
        :rtype: tuple

        """

        # Create a polygon of the envelope of the first image
        xmin, xmax, ymin, ymax = self.extent
        wkt = "POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))" \
            %(xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin,xmin,ymin)
        poly1 = ogr.CreateGeometryFromWkt(wkt)

        # Create a polygon of the envelope of the second image
        img = SingleBandRaster(filename, load_data=False)
        xmin, xmax, ymin, ymax = img.extent
        wkt = "POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))" \
            %(xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin,xmin,ymin)
        poly2 = ogr.CreateGeometryFromWkt(wkt)

        # Compute intersection envelope
        intersect = poly1.Intersection(poly2)
        extent = intersect.GetEnvelope()

        # check that intersection is not void
        if intersect.GetArea()==0:
            print('Warning: Intersection is void')
            return 0
        else:
            return extent

def has_geos():
    pnt1 = ogr.CreateGeometryFromWkt('POINT(10 20)')
    pnt2 = ogr.CreateGeometryFromWkt('POINT(30 20)')
    ogrex = ogr.GetUseExceptions()
    ogr.DontUseExceptions()
    hasgeos = pnt1.Union(pnt2) is not None
    if ogrex:
        ogr.UseExceptions()
    return hasgeos

def copy_features(src_layer, dst_layer, fix_geometry, simplify_geometry, start_index, total):
    index = 0
    batch_size = 200
    index_batch = 0
    for feat in src_layer:
        if index < start_index:
            index = index + 1
            continue

        try:
            geom = shapely.wkt.loads(feat.GetGeometryRef().ExportToWkt())
        except Exception as e:
            print('Error({0}), skipping geometry.'.format(e))
            continue

        if fix_geometry and not geom.is_valid:
            geom = geom.buffer(0.0)

        if simplify_geometry:
            geom = geom.simplify(0.004)

        f = ogr.Feature(dst_layer.GetLayerDefn())

        # set field values
        for i in range(feat.GetFieldCount()):
            fd = feat.GetFieldDefnRef(i)
            f.SetField(fd.GetName(), feat.GetField(fd.GetName()))

        # set geometry    
        f.SetGeometry(ogr.CreateGeometryFromWkt(geom.to_wkt()))

        if index_batch == 0:
            dst_layer.StartTransaction()

        # create feature
        feature = dst_layer.CreateFeature(f)

        f.Destroy()

        index_batch = index_batch + 1

        if index_batch >= batch_size or index == total - 1:
            dst_layer.CommitTransaction()
            count = dst_layer.GetFeatureCount() # update number of inserted features
            print('Inserted {0} of {1} features ({2:.2f}%)'.format(count, total, 100. * float(count) / total))

            index_batch = 0

            if index == total - 1:
                break

        index = index + 1

def shp2array(shp_fn, r_ds=None, res=None, extent=None, t_srs=None):
    """Rasterize input shapefile to match existing raster Dataset (or specified res/extent/t_srs)
    """
    shp_ds = ogr.Open(shp_fn)
    lyr = shp_ds.GetLayer()
    #This returns xmin, ymin, xmax, ymax
    shp_extent = lyr_extent(lyr)
    shp_srs = lyr.GetSpatialRef()
    # dst_dt = gdal.GDT_Byte
    ndv = 0
    if r_ds is not None:
        r_extent = ds_extent(r_ds)
        res = get_res(r_ds, square=True)[0] 
        if extent is None:
            extent = r_extent
        r_srs = get_ds_srs(r_ds)
        r_geom = ds_geom(r_ds)
        # dst_ns = r_ds.RasterXSize
        # dst_nl = r_ds.RasterYSize
        #Convert raster extent to shp_srs
        cT = osr.CoordinateTransformation(r_srs, shp_srs)
        r_geom_reproj = geom_dup(r_geom)
        r_geom_reproj.Transform(cT)
        r_geom_reproj.AssignSpatialReference(t_srs)
        lyr.SetSpatialFilter(r_geom_reproj)
        #lyr.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt))
    else:
        #TODO: clean this up
        if res is None:
            sys.exit("Must specify input res")
        if extent is None:
            print("Using input shp extent")
            extent = shp_extent
    if t_srs is None:
        t_srs = r_srs
    if not shp_srs.IsSame(t_srs):
        print("Input shp srs: %s" % shp_srs.ExportToProj4())
        print("Specified output srs: %s" % t_srs.ExportToProj4())
        out_ds = lyr_proj(lyr, t_srs)
        outlyr = out_ds.GetLayer()
    else:
        outlyr = lyr
    #outlyr.SetSpatialFilter(r_geom)
    m_ds = mem_ds(res, extent, srs=t_srs, dtype=gdal.GDT_Byte)
    b = m_ds.GetRasterBand(1)
    b.SetNoDataValue(ndv)
    gdal.RasterizeLayer(m_ds, [1], outlyr, burn_values=[1])
    a = b.ReadAsArray()
    a = ~(a.astype('Bool'))
    return a