Python osgeo.ogr.CreateGeometryFromWkt() Examples

The following are 15 code examples of osgeo.ogr.CreateGeometryFromWkt(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module osgeo.ogr , or try the search function .
Example #1
Source Project: momepy   Author: martinfleis   File: elements.py    License: MIT License 7 votes vote down vote up
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 
Example #2
Source Project: cloudless   Author: BradNeuberg   File: download_planetlabs.py    License: Apache License 2.0 6 votes vote down vote up
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) 
Example #3
Source Project: unmixing   Author: arthur-e   File: lsma.py    License: MIT License 6 votes vote down vote up
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() 
Example #4
Source Project: neon   Author: NervanaSystems   File: spacenet_utils.py    License: Apache License 2.0 6 votes vote down vote up
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 
Example #5
Source Project: buzzard   Author: airware   File: _gdal_conv.py    License: Apache License 2.0 5 votes vote down vote up
def ogr_of_shapely(geom):
    return ogr.CreateGeometryFromWkt(geom.wkt) 
Example #6
Source Project: cloudless   Author: BradNeuberg   File: download_planetlabs.py    License: Apache License 2.0 5 votes vote down vote up
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 
Example #7
Source Project: pyeo   Author: clcr   File: coordinate_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
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 
Example #8
Source Project: cadasta-platform   Author: Cadasta   File: shape.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
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() 
Example #9
Source Project: esa_sentinel   Author: jonas-eberle   File: sentinel_api.py    License: MIT License 5 votes vote down vote up
def multipolygon2list(wkt):
        geom = ogr.CreateGeometryFromWkt(wkt)
        if geom.GetGeometryName() == 'MULTIPOLYGON':
            return [x.ExportToWkt() for x in geom]
        else:
            return [geom.ExportToWkt()] 
Example #10
Source Project: wradlib   Author: wradlib   File: util.py    License: MIT License 5 votes vote down vote up
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 
Example #11
Source Project: wradlib   Author: wradlib   File: test_georef.py    License: MIT License 5 votes vote down vote up
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)) 
Example #12
Source Project: pyTSEB   Author: hectornieto   File: PyTSEB.py    License: GNU General Public License v3.0 4 votes vote down vote up
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 
Example #13
Source Project: geoio   Author: DigitalGlobe   File: base.py    License: MIT License 4 votes vote down vote up
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.") 
Example #14
Source Project: gips   Author: gipit   File: core.py    License: GNU General Public License v2.0 4 votes vote down vote up
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 
Example #15
Source Project: RHEAS   Author: nasa   File: analysis.py    License: MIT License 4 votes vote down vote up
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()