Python osgeo.ogr.CreateGeometryFromWkt() Examples

def _densify(self, geom, segment):
        """
        Returns densified geoemtry with segments no longer than `segment`.
        """
        # temporary solution for readthedocs fail. - cannot mock osgeo
        try:
            from osgeo import ogr
        except ModuleNotFoundError:
            import warnings

            warnings.warn("OGR (GDAL) is required.")

        poly = geom
        wkt = geom.wkt  # shapely Polygon to wkt
        geom = ogr.CreateGeometryFromWkt(wkt)  # create ogr geometry
        geom.Segmentize(segment)  # densify geometry by 2 metres
        geom.CloseRings()  # fix for GDAL 2.4.1 bug
        wkt2 = geom.ExportToWkt()  # ogr geometry to wkt
        try:
            new = loads(wkt2)  # wkt to shapely Polygon
            return new
        except Exception:
            return poly 

def buffer_bbox(geom, buff):
    """
    Buffers the geom by buff and then calculates the bounding box.
    Returns a Geometry of the bounding box
    """
    b = geom.Buffer(buff)
    lng1, lng2, lat1, lat2 = b.GetEnvelope()
    wkt = """POLYGON((
        %s %s,
        %s %s,
        %s %s,
        %s %s,
        %s %s
    ))""" % (lng1, lat1, lng1, lat2, lng2, lat2, lng2, lat1, lng1, lat1)
    wkt = wkt.replace('\n', '')
    return ogr.CreateGeometryFromWkt(wkt) 

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 get_bounding_boxes(img_file, annot_file):

    srcRaster = gdal.Open(img_file)
    targetSR = osr.SpatialReference()
    targetSR.ImportFromWkt(srcRaster.GetProjectionRef())
    geomTransform = srcRaster.GetGeoTransform()

    dataSource = ogr.Open(annot_file, 0)
    layer = dataSource.GetLayer()

    building_id = 0
    buildinglist = []

    for feature in layer:
        geom = feature.GetGeometryRef()
        geom_wkt_list = geoPolygonToPixelPolygonWKT(geom, img_file, targetSR, geomTransform)

        for geom_wkt in geom_wkt_list:
            building_id += 1
            buildinglist.append(ogr.CreateGeometryFromWkt(geom_wkt[0]).GetEnvelope())

    return buildinglist 

def ogr_of_shapely(geom):
    return ogr.CreateGeometryFromWkt(geom.wkt) 

def download(lat, lng, buff_meters, download_dir='/tmp', image_type='planetlabs'):
    """
    Given a latitude, longitude, and a number of meters to buffer by,
    get all imagery that intersects the bounding box of the buffer and
    download it to the specified directory.  Return the names of the
    downloaded files.
    """
    pt = ogr.CreateGeometryFromWkt('POINT(%s %s)' % (lng, lat))
    pt = reproject(pt, 4326, 2163)  # from WGS84 to National Atlas
    buff = buffer_bbox(pt, buff_meters)
    buff = reproject(buff, 2163, 4326)

    if image_type == 'planetlabs':
        scenes_url = "https://api.planet.com/v0/scenes/ortho/"
    elif image_type == 'rapideye':
        scenes_url = "https://api.planet.com/v0/scenes/rapideye/"

    # Download the initial scenes URL and also any paged results that come after that.
    downloaded_scenes = []
    next_url = scenes_url
    params = {"intersects": buff.ExportToWkt()}
    while next_url != None:
        next_url = download_results(next_url, params, downloaded_scenes, download_dir, image_type)
        print "\nWorking with next page of results: %s" % next_url

    return downloaded_scenes 

def point_to_pixel_coordinates(raster, point, oob_fail=False):
    """
    Returns a tuple (x_pixel, y_pixel) in a georaster raster corresponding to the geographic point in a projection.
     Assumes raster is north-up non rotated.

    Parameters
    ----------
    raster
        A gdal raster object
    point
        One of:
            A well-known text string of a single point
            An iterable of the form (x,y)
            An ogr.Geometry object containing a single point
    Returns
    -------
    A tuple of (x_pixel, y_pixel), containing the indicies of the point in the raster.

    Notes
    -----
    The equation is a rearrangement of the section on affinine geotransform in http://www.gdal.org/gdal_datamodel.html

    """
    if isinstance(point, str):
        point = ogr.CreateGeometryFromWkt(point)
        x_geo = point.GetX()
        y_geo = point.GetY()
    if isinstance(point, list) or isinstance(point, tuple):  # There is a more pythonic way to do this
        x_geo = point[0]
        y_geo = point[1]
    if isinstance(point, ogr.Geometry):
        x_geo = point.GetX()
        y_geo = point.GetY()
    gt = raster.GetGeoTransform()
    x_pixel = int(np.floor((x_geo - floor_to_resolution(gt[0], gt[1]))/gt[1]))
    y_pixel = int(np.floor((y_geo - floor_to_resolution(gt[3], gt[5]*-1))/gt[5]))  # y resolution is -ve
    return x_pixel, y_pixel 

def write_shp_layer(self, loc_data):
        if loc_data['geometry.wkt'] == '':
            return
        geom = ogr.CreateGeometryFromWkt(loc_data['geometry.wkt'])
        layer_type = geom.GetGeometryName().lower()
        layer = self.shp_layers.get(layer_type, None)
        if layer is None:
            layer = self.create_shp_layer(layer_type)
            self.shp_layers[layer_type] = layer
        if layer:
            feature = ogr.Feature(layer.GetLayerDefn())
            feature.SetGeometry(geom)
            feature.SetField('id', loc_data['id'])
            layer.CreateFeature(feature)
            feature.Destroy() 

def multipolygon2list(wkt):
        geom = ogr.CreateGeometryFromWkt(wkt)
        if geom.GetGeometryName() == 'MULTIPOLYGON':
            return [x.ExportToWkt() for x in geom]
        else:
            return [geom.ExportToWkt()] 

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

def test_get_vector_points_warning(self):
        point_wkt = "POINT (1198054.34 648493.09)"
        point = ogr.CreateGeometryFromWkt(point_wkt)
        with pytest.warns(UserWarning):
            list(georef.get_vector_points(point)) 

def _get_subset(self, roi_shape, raster_proj_wkt, raster_geo_transform):

        # Find extent of ROI in roiShape projection
        roi = ogr.Open(roi_shape)
        roi_layer = roi.GetLayer()
        roi_extent = roi_layer.GetExtent()

        # Convert the extent to raster projection
        roi_proj = roi_layer.GetSpatialRef()
        raster_proj = osr.SpatialReference()
        raster_proj.ImportFromWkt(raster_proj_wkt)
        transform = osr.CoordinateTransformation(roi_proj, raster_proj)
        point_UL = ogr.CreateGeometryFromWkt("POINT ({} {})"
                                             .format(min(roi_extent[0], roi_extent[1]),
                                                     max(roi_extent[2], roi_extent[3])))
        point_UL.Transform(transform)
        point_UL = point_UL.GetPoint()
        point_LR = ogr.CreateGeometryFromWkt("POINT ({} {})"
                                             .format(max(roi_extent[0], roi_extent[1]),
                                                     min(roi_extent[2], roi_extent[3])))
        point_LR.Transform(transform)
        point_LR = point_LR.GetPoint()

        # Get pixel location of this extent
        ulX = raster_geo_transform[0]
        ulY = raster_geo_transform[3]
        pixel_size = raster_geo_transform[1]
        pixel_UL = [max(int(math.floor((ulY - point_UL[1]) / pixel_size)), 0),
                    max(int(math.floor((point_UL[0] - ulX) / pixel_size)), 0)]
        pixel_LR = [int(round((ulY - point_LR[1]) / pixel_size)),
                    int(round((point_LR[0] - ulX) / pixel_size))]

        # Get projected extent
        point_proj_UL = (ulX + pixel_UL[1] * pixel_size, ulY - pixel_UL[0] * pixel_size)

        # Convert to xoff, yoff, xcount, ycount as required by GDAL ReadAsArray()
        subset_pix = [pixel_UL[1], pixel_UL[0],
                      pixel_LR[1] - pixel_UL[1], pixel_LR[0] - pixel_UL[0]]

        # Get the geo transform of the subset
        subset_geo_transform = [point_proj_UL[0], pixel_size, raster_geo_transform[2],
                                point_proj_UL[1], raster_geo_transform[4], -pixel_size]

        return subset_pix, subset_geo_transform 

def _instantiate_geom(self,g):
        """Attempt to convert the geometry pass in to and ogr Geometry
        object.  Currently implements the base ogr.CreateGeometryFrom*
        methods and will reform fiona geometry dictionaries into a format
        that ogr.CreateGeometryFromJson will correctly handle.
        """

        if isinstance(g,ogr.Geometry):
            # If the input geometry is already an ogr object, create a copy
            # of it.  This is requred because of a subtle issue that causes
            # gdal to crash if the input geom is used elsewhere.  The working
            # theory is that the geometry is freed when going out of a scope
            # while it is needed in the upper level loop.  In this code, the
            # problem comes about between self.iter_vector and self.get_data
            # with mask=True.
            return ogr.CreateGeometryFromJson(str(g.ExportToJson()))

        # Handle straight ogr GML
        try:
            return ogr.CreateGeometryFromGML(g)
        except:
            pass
        # Handle straight ogr Wkb
        try:
            return ogr.CreateGeometryFromWkb(g)
        except:
            pass
        # Handle straight ogr Json
        try:
            return ogr.CreateGeometryFromJson(g)
        except:
            pass
        # Handle straight ogr Wkt
        try:
            return ogr.CreateGeometryFromWkt(g)
        except:
            pass
        # Handle fiona Json geometry format
        try:
            gjson = str(g).replace('(','[').replace(')',']')
            return ogr.CreateGeometryFromJson(gjson)
        except:
            pass

        raise ValueError("A geometry object was not able to be created from " \
                         "the value passed in.") 

def vector2tiles(cls, vector, pcov=0.0, ptile=0.0, tilelist=None):
        """ Return matching tiles and coverage % for provided vector """
        from osgeo import ogr, osr

        # open tiles vector
        v = open_vector(cls.get_setting('tiles'))
        shp = ogr.Open(v.Filename())
        if v.LayerName() == '':
            layer = shp.GetLayer(0)
        else:
            layer = shp.GetLayer(v.LayerName())

        # create and warp site geometry
        ogrgeom = ogr.CreateGeometryFromWkt(vector.WKT())
        srs = osr.SpatialReference(vector.Projection())
        trans = osr.CoordinateTransformation(srs, layer.GetSpatialRef())
        ogrgeom.Transform(trans)
        # convert to shapely
        geom = loads(ogrgeom.ExportToWkt())

        # find overlapping tiles
        tiles = {}
        layer.SetSpatialFilter(ogrgeom)
        layer.ResetReading()
        feat = layer.GetNextFeature()
        while feat is not None:
            tgeom = loads(feat.GetGeometryRef().ExportToWkt())
            if tgeom.intersects(geom):
                area = geom.intersection(tgeom).area
                if area != 0:
                    tile = cls.feature2tile(feat)
                    tiles[tile] = (area / geom.area, area / tgeom.area)
            feat = layer.GetNextFeature()

        # remove any tiles not in tilelist or that do not meet thresholds for % cover
        remove_tiles = []
        if tilelist is None:
            tilelist = tiles.keys()
        for t in tiles:
            if (tiles[t][0] < (pcov / 100.0)) or (tiles[t][1] < (ptile / 100.0)) or t not in tilelist:
                remove_tiles.append(t)
        for t in remove_tiles:
            tiles.pop(t, None)
        return tiles 

def cropYield(shapefile, name, startdate="", enddate="", crop="maize", dbname="rheas"):
    """Extract crop yield from a specified simulation *name* for dates ranging
    from *startdate* to *enddate*, and saves them a *shapefile*."""
    logging.basicConfig(level=logging.INFO, format='%(message)s')
    log = logging.getLogger(__name__)
    db = dbio.connect(dbname)
    cur = db.cursor()
    datesql = ""
    if len(startdate) > 0:
        try:
            sdt = datetime.strptime(startdate, "%Y-%m-%d")
            datesql = "and fdate>=date'{0}'".format(sdt.strftime("%Y-%m-%d"))
        except ValueError:
            log.warning("Start date is invalid and will be ignored.")
    if len(enddate) > 0:
        try:
            edt = datetime.strptime(enddate, "%Y-%m-%d")
            datesql += "and fdate<=date'{0}'".format(edt.strftime("%y-%m-%d"))
        except ValueError:
            log.warning("End date is invalid and will be ignored.")
    fsql = "with f as (select gid,geom,gwad,ensemble,fdate from (select gid,geom,gwad,ensemble,fdate,row_number() over (partition by gid,ensemble order by gwad desc) as rn from {0}.dssat) gwadtable where rn=1 {1})".format(name, datesql)
    sql = "{0} select gid,st_astext(geom),max(gwad) as max_yield,avg(gwad) as avg_yield,stddev(gwad) as std_yield,max(fdate) as fdate from f group by gid,geom".format(fsql)
    cur.execute(sql)
    if bool(cur.rowcount):
        results = cur.fetchall()
        drv = ogr.GetDriverByName("ESRI Shapefile")
        ds = drv.CreateDataSource(shapefile)
        lyr = ds.CreateLayer("yield", geom_type=ogr.wkbMultiPolygon)
        lyr.CreateField(ogr.FieldDefn("gid", ogr.OFTInteger))
        lyr.CreateField(ogr.FieldDefn("average", ogr.OFTReal))
        lyr.CreateField(ogr.FieldDefn("maximum", ogr.OFTReal))
        lyr.CreateField(ogr.FieldDefn("minimum", ogr.OFTReal))
        for row in results:
            feat = ogr.Feature(lyr.GetLayerDefn())
            feat.SetField("gid", row[0])
            feat.SetField("maximum", row[2])
            feat.SetField("average", row[3])
            feat.SetField("minimum", row[4])
            feat.SetGeometry(ogr.CreateGeometryFromWkt(row[1]))
            lyr.CreateFeature(feat)
            feat.Destroy()
        ds.Destroy()