Python osgeo.ogr.Open() Examples

The following are 30 code examples of osgeo.ogr.Open(). 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: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 7 votes vote down vote up
def clip_raster_to_intersection(raster_to_clip_path, extent_raster_path, out_raster_path, is_landsat=False):
    """
    Clips one raster to the extent proivded by the other raster, and saves the result at out_raster_path.
    Assumes both raster_to_clip and extent_raster are in the same projection.
    Parameters
    ----------
    raster_to_clip_path
        The location of the raster to be clipped.
    extent_raster_path
        The location of the raster that will provide the extent to clip to
    out_raster_path
        A location for the finished raster
    """

    with TemporaryDirectory() as td:
        temp_aoi_path = os.path.join(td, "temp_clip.shp")
        get_extent_as_shp(extent_raster_path, temp_aoi_path)
        ext_ras = gdal.Open(extent_raster_path)
        proj = osr.SpatialReference(wkt=ext_ras.GetProjection())
        srs_id = int(proj.GetAttrValue('AUTHORITY', 1))
        clip_raster(raster_to_clip_path, temp_aoi_path, out_raster_path, srs_id, flip_x_y = is_landsat) 
Example #2
Source Project: TF-SegNet   Author: mathildor   File: import_data.py    License: MIT License 7 votes vote down vote up
def importAr(to_cur, filename):
    print("Importing:", filename)
    dataSource = ogr.Open(filename)
    dataLayer = dataSource.GetLayer(0)

    for feature in dataLayer:
        geom = feature.GetGeometryRef().ExportToWkt()
        id = feature.GetField("poly_id")
        objtype = feature.GetField("objtype")
        artype = feature.GetField("artype")
        arskogbon = feature.GetField("arskogbon")
        artreslag = feature.GetField("artreslag")
        argrunnf = feature.GetField("argrunnf")

        to_tuple = ( id, objtype, artype, arskogbon, artreslag, argrunnf, geom)
        to_cur.execute("""INSERT into ar_bygg (id, objtype, artype, arskogbon, artreslag, argrunnf, geom)
                        SELECT %s, %s, %s, %s, %s, %s, ST_GeometryFromText(%s);""",
                   to_tuple)

    to_conn.commit()
    dataSource.Destroy() 
Example #3
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def strip_bands(in_raster_path, out_raster_path, bands_to_strip):
    in_raster = gdal.Open(in_raster_path)
    out_raster_band_count = in_raster.RasterCount-len(bands_to_strip)
    out_raster = create_matching_dataset(in_raster, out_raster_path, bands=out_raster_band_count)
    out_raster_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    in_raster_array = in_raster.GetVirtualMemArray()

    bands_to_copy = [band for band in range(in_raster_array.shape[0]) if band not in bands_to_strip]

    out_raster_array[...] = in_raster_array[bands_to_copy, :,:]

    out_raster_array = None
    in_raster_array = None
    out_raster = None
    in_raster = None

    return out_raster_path 
Example #4
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def flatten_probability_image(prob_image, out_path):
    """
    Takes a probability output from classify_image and flattens it into a single layer containing only the maximum
    value from each pixel.

    Parameters
    ----------
    prob_image
        The path to a probability image.
    out_path
        The place to save the flattened image.

    """
    prob_raster = gdal.Open(prob_image)
    out_raster = create_matching_dataset(prob_raster, out_path, bands=1)
    prob_array = prob_raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[:, :] = prob_array.max(axis=0)
    out_array = None
    prob_array = None
    out_raster = None
    prob_raster = None 
Example #5
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def raster_to_array(rst_pth):
    """Reads in a raster file and returns a N-dimensional array.

    Parameters
    ----------
    rst_pth
        Path to input raster.
    Returns
    -------
    As N-dimensional array.
    """
    log = logging.getLogger(__name__)
    in_ds = gdal.Open(rst_pth)
    out_array = in_ds.ReadAsArray()

    return out_array 
Example #6
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def calc_ndvi(raster_path, output_path):
    raster = gdal.Open(raster_path)
    out_raster = create_matching_dataset(raster, output_path, datatype=gdal.GDT_Float32)
    array = raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    R = array[2, ...]
    I = array[3, ...]
    out_array[...] = (R-I)/(R+I)

    out_array[...] = np.where(out_array == -2147483648, 0, out_array)

    R = None
    I = None
    array = None
    out_array = None
    raster = None
    out_raster = None 
Example #7
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def apply_image_function(in_paths, out_path, function, out_datatype = gdal.GDT_Int32):
    """Applies a pixel-wise function across every image. Assumes each image is exactly contiguous and, for now,
    single-banded. function() should take a list of values and return a single value."""
    rasters = [gdal.Open(in_path) for in_path in in_paths]
    raster_arrays = [raster.GetVirtualMemArray() for raster in rasters]
    in_array = np.stack(raster_arrays, axis=0)

    out_raster = create_matching_dataset(rasters[0], out_path=out_path, datatype=out_datatype)
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[...] = np.apply_along_axis(function, 0, in_array)

    # Deallocating. Not taking any chances here.
    out_array = None
    out_raster = None
    in_array = None
    for raster_array, raster in zip(raster_arrays, rasters):
        raster_array = None
        raster = None 
Example #8
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def create_mask_from_model(image_path, model_path, model_clear=0, num_chunks=10, buffer_size=0):
    """Returns a multiplicative mask (0 for cloud, shadow or haze, 1 for clear) built from the model at model_path."""
    from pyeo.classification import classify_image  # Deferred import to deal with circular reference
    with TemporaryDirectory() as td:
        log = logging.getLogger(__name__)
        log.info("Building cloud mask for {} with model {}".format(image_path, model_path))
        temp_mask_path = os.path.join(td, "cat_mask.tif")
        classify_image(image_path, model_path, temp_mask_path, num_chunks=num_chunks)
        temp_mask = gdal.Open(temp_mask_path, gdal.GA_Update)
        temp_mask_array = temp_mask.GetVirtualMemArray()
        mask_path = get_mask_path(image_path)
        mask = create_matching_dataset(temp_mask, mask_path, datatype=gdal.GDT_Byte)
        mask_array = mask.GetVirtualMemArray(eAccess=gdal.GF_Write)
        mask_array[:, :] = np.where(temp_mask_array != model_clear, 0, 1)
        temp_mask_array = None
        mask_array = None
        temp_mask = None
        mask = None
        if buffer_size:
            buffer_mask_in_place(mask_path, buffer_size)
        log.info("Cloud mask for {} saved in {}".format(image_path, mask_path))
        return mask_path 
Example #9
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def create_mask_from_class_map(class_map_path, out_path, classes_of_interest, buffer_size=0, out_resolution=None):
    """Creates a mask from a classification mask: 1 for each pixel containing one of classes_of_interest, otherwise 0"""
    # TODO: pull this out of the above function
    class_image = gdal.Open(class_map_path)
    class_array = class_image.GetVirtualMemArray()
    mask_array = np.isin(class_array, classes_of_interest)
    out_mask = create_matching_dataset(class_image, out_path, datatype=gdal.GDT_Byte)
    out_array = out_mask.GetVirtualMemArray(eAccess=gdal.GA_Update)
    np.copyto(out_array, mask_array)
    class_array = None
    class_image = None
    out_array = None
    out_mask = None
    if out_resolution:
        resample_image_in_place(out_path, out_resolution)
    if buffer_size:
        buffer_mask_in_place(out_path, buffer_size)
    return out_path 
Example #10
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 6 votes vote down vote up
def create_mask_from_fmask(in_l1_dir, out_path):
    log = logging.getLogger(__name__)
    log.info("Creating fmask for {}".format(in_l1_dir))
    with TemporaryDirectory() as td:
        temp_fmask_path = os.path.join(td, "fmask.tif")
        apply_fmask(in_l1_dir, temp_fmask_path)
        fmask_image = gdal.Open(temp_fmask_path)
        fmask_array = fmask_image.GetVirtualMemArray()
        out_image = create_matching_dataset(fmask_image, out_path, datatype=gdal.GDT_Byte)
        out_array = out_image.GetVirtualMemArray(eAccess=gdal.GA_Update)
        log.info("fmask created, converting to binary cloud/shadow mask")
        out_array[:,:] = np.isin(fmask_array, (2, 3, 4), invert=True)
        out_array = None
        out_image = None
        fmask_array = None
        fmask_image = None
        resample_image_in_place(out_path, 10) 
Example #11
Source Project: CityEnergyAnalyst   Author: architecture-building-systems   File: terrain_helper.py    License: MIT License 6 votes vote down vote up
def calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, x_max, x_min, y_max,
                                        y_min):
    # local variables:
    temp_shapefile = locator.get_temporary_file("terrain.shp")
    cols = int((x_max - x_min) / grid_size)
    rows = int((y_max - y_min) / grid_size)
    shapes = Polygon([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max], [x_min, y_min]])
    geodataframe = Gdf(index=[0], crs=crs, geometry=[shapes])
    geodataframe.to_file(temp_shapefile)
    # 1) opening the shapefile
    source_ds = ogr.Open(temp_shapefile)
    source_layer = source_ds.GetLayer()
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Float32)  ##COMMENT 2
    target_ds.SetGeoTransform((x_min, grid_size, 0, y_max, 0, -grid_size))  ##COMMENT 3
    # 5) Adding a spatial reference ##COMMENT 4
    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromProj4(crs)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(-9999)  ##COMMENT 5
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[elevation])  ##COMMENT 6
    target_ds = None  # closing the file 
Example #12
Source Project: osgeopy-code   Author: cgarrard   File: listing4_4.py    License: MIT License 6 votes vote down vote up
def add_markers(fmap, json_fn):
    ds = ogr.Open(json_fn)
    lyr = ds.GetLayer()
    for row in lyr:
        geom = row.geometry()
        color = colors[row.GetField('status')]
        fmap.circle_marker([geom.GetY(), geom.GetX()],
                           line_color=color,
                           fill_color=color,
                           radius=5000,
                           popup=get_popup(row.items()))

# Because this version of make_map was defined after the version in
# listing4_3 (since it was already imported and evaluated), then this is
# the version that will be used. Unlike in the book text, the first two
# lines of the function here have been changed to reference get_state_geom,
# save_state_gauges, get_bbox, and get_center from the listing4_3 module. 
Example #13
Source Project: osgeopy-code   Author: cgarrard   File: listing4_3.py    License: MIT License 6 votes vote down vote up
def save_state_gauges(out_fn, bbox=None):
    """Save stream gauge data to a geojson file."""
    url = 'http://gis.srh.noaa.gov/arcgis/services/ahps_gauges/' + \
          'MapServer/WFSServer'
    parms = {
        'version': '1.1.0',
        'typeNames': 'ahps_gauges:Observed_River_Stages',
        'srsName': 'urn:ogc:def:crs:EPSG:6.9:4326',
    }
    if bbox:
        parms['bbox'] = bbox
    try:
        request = 'WFS:{0}?{1}'.format(url, urllib.urlencode(parms))
    except:
        request = 'WFS:{0}?{1}'.format(url, urllib.parse.urlencode(parms))
    wfs_ds = ogr.Open(request)
    if wfs_ds is None:
        raise RuntimeError('Could not open WFS.')
    wfs_lyr = wfs_ds.GetLayer(0)

    driver = ogr.GetDriverByName('GeoJSON')
    if os.path.exists(out_fn):
        driver.DeleteDataSource(out_fn)
    json_ds = driver.CreateDataSource(out_fn)
    json_ds.CopyLayer(wfs_lyr, '') 
Example #14
Source Project: stdm   Author: gltn   File: gps_tool_data_source_utils.py    License: GNU General Public License v2.0 6 votes vote down vote up
def open_gpx_file(gpx_file):
    """
    Opens input GPX file
    :param gpx_file: Input GPX file
    :return: File object
    :rtype: Object
    """
    gpx_data_source = ogr.Open(gpx_file)
    if not gpx_data_source:
        raise Exception(
            "File {0} could not be accessed. "
            "It could be in use by another application".format(
                os.path.basename(gpx_file)
            )
        )
    return gpx_data_source 
Example #15
Source Project: stdm   Author: gltn   File: gps_tool_data_source_utils.py    License: GNU General Public License v2.0 6 votes vote down vote up
def get_feature_layer(gpx_data_source, feature_type):
    """
    Gets valid feature layer from the GPX file based on
    the selected feature type
    :param gpx_data_source: Open file object
    :param feature_type: User feature type input
    :return  gpx_layer: GPX file feature layer
    :return feature_count: Number of features
    :rtype gpx_layer: Layer object
    :rtype feature_count: Integer
    """
    if feature_type >= 0:
        if FEATURE_TYPES[feature_type] == 'track' or \
                        FEATURE_TYPES[feature_type] == 'route':
            feature_type = '{}_points'.format(FEATURE_TYPES[feature_type])
        else:
            feature_type = '{}s'.format(FEATURE_TYPES[feature_type])
        gpx_layer = gpx_data_source.GetLayerByName(feature_type)
        feature_count = gpx_layer.GetFeatureCount()
        return None if feature_count == 0 else gpx_layer, feature_count 
Example #16
Source Project: wa   Author: wateraccounting   File: raster_conversions.py    License: Apache License 2.0 6 votes vote down vote up
def Open_array_info(filename=''):
    """
    Opening a tiff info, for example size of array, projection and transform matrix.

    Keyword Arguments:
    filename -- 'C:/file/to/path/file.tif' or a gdal file (gdal.Open(filename))
        string that defines the input tiff file or gdal file

    """
    f = gdal.Open(r"%s" %filename)
    if f is None:
        print '%s does not exists' %filename
    else:
        geo_out = f.GetGeoTransform()
        proj = f.GetProjection()
        size_X = f.RasterXSize
        size_Y = f.RasterYSize
        f = None
    return(geo_out, proj, size_X, size_Y) 
Example #17
Source Project: wa   Author: wateraccounting   File: raster_conversions.py    License: Apache License 2.0 6 votes vote down vote up
def Open_tiff_array(filename='', band=''):
    """
    Opening a tiff array.

    Keyword Arguments:
    filename -- 'C:/file/to/path/file.tif' or a gdal file (gdal.Open(filename))
        string that defines the input tiff file or gdal file
    band -- integer
        Defines the band of the tiff that must be opened.
    """
    f = gdal.Open(filename)
    if f is None:
        print '%s does not exists' %filename
    else:
        if band is '':
            band = 1
        Data = f.GetRasterBand(band).ReadAsArray()
    return(Data) 
Example #18
Source Project: DsgTools   Author: dsgoficial   File: inventoryThread.py    License: GNU General Public License v2.0 6 votes vote down vote up
def copy(self, destinationFolder):
        """
        Copy inventoried files considering the dataset
        destinationFolder: copy destination folder
        """
        for fileName in self.files:
            # adjusting the separators according to the OS
            fileName = fileName.replace('/', os.sep)
            file = fileName.split(os.sep)[-1]
            newFileName = os.path.join(destinationFolder, file)
            newFileName = newFileName.replace('/', os.sep)

            try:
                gdalSrc = gdal.Open(fileName)
                ogrSrc = ogr.Open(fileName)
                if ogrSrc:
                    self.copyOGRDataSource(ogrSrc, newFileName)
                elif gdalSrc:
                    self.copyGDALDataSource(gdalSrc, newFileName)
            except Exception as e:
                QgsMessageLog.logMessage(self.messenger.getCopyErrorMessage()+'\n'+':'.join(e.args), "DSGTools Plugin", QgsMessageLog.INFO)
                return (0, self.messenger.getCopyErrorMessage()+'\n'+':'.join(e.args))
        
        QgsMessageLog.logMessage(self.messenger.getSuccessInventoryAndCopyMessage(), "DSGTools Plugin", QgsMessageLog.INFO)
        return (1, self.messenger.getSuccessInventoryAndCopyMessage()) 
Example #19
Source Project: gips   Author: gipit   File: core.py    License: GNU General Public License v2.0 6 votes vote down vote up
def datafiles(self):
        """ Get list of readable datafiles from asset (multiple filenames if tar or hdf file) """
        path = os.path.dirname(self.filename)
        indexfile = os.path.join(path, self.filename + '.index')
        if os.path.exists(indexfile):
            datafiles = File2List(indexfile)
            if len(datafiles) > 0:
                return datafiles
        try:
            if tarfile.is_tarfile(self.filename):
                tfile = tarfile.open(self.filename)
                tfile = tarfile.open(self.filename)
                datafiles = tfile.getnames()
            else:
                # Try subdatasets
                fh = gdal.Open(self.filename)
                sds = fh.GetSubDatasets()
                datafiles = [s[0] for s in sds]
            if len(datafiles) > 0:
                List2File(datafiles, indexfile)
                return datafiles
            else:
                return [self.filename]
        except:
            raise Exception('Problem accessing asset(s) in %s' % self.filename) 
Example #20
Source Project: qgisSpaceSyntaxToolkit   Author: SpaceGroupUCL   File: test_shp.py    License: GNU General Public License v3.0 6 votes vote down vote up
def test_attributeexport(self):
        def testattributes(lyr, graph):
            feature = lyr.GetNextFeature()
            while feature:
                coords = []
                ref = feature.GetGeometryRef()
                last = ref.GetPointCount() - 1
                edge_nodes = (ref.GetPoint_2D(0), ref.GetPoint_2D(last))
                name = feature.GetFieldAsString('Name')
                assert_equal(graph.get_edge_data(*edge_nodes)['Name'], name)
                feature = lyr.GetNextFeature()

        tpath = os.path.join(tempfile.gettempdir(), 'shpdir')

        G = nx.read_shp(self.shppath)
        nx.write_shp(G, tpath)
        shpdir = ogr.Open(tpath)
        edges = shpdir.GetLayerByName("edges")
        testattributes(edges, G) 
Example #21
Source Project: qgisSpaceSyntaxToolkit   Author: SpaceGroupUCL   File: test_shp.py    License: GNU General Public License v3.0 6 votes vote down vote up
def test_wkt_export(self):
        G = nx.DiGraph()
        tpath = os.path.join(tempfile.gettempdir(), 'shpdir')
        points = (
            "POINT (0.9 0.9)",
            "POINT (4 2)"
        )
        line = (
            "LINESTRING (0.9 0.9,4 2)",
        )
        G.add_node(1, Wkt=points[0])
        G.add_node(2, Wkt=points[1])
        G.add_edge(1, 2, Wkt=line[0])
        try:
            nx.write_shp(G, tpath)
        except Exception as e:
            assert False, e
        shpdir = ogr.Open(tpath)
        self.checkgeom(shpdir.GetLayerByName("nodes"), points)
        self.checkgeom(shpdir.GetLayerByName("edges"), line) 
Example #22
Source Project: RHEAS   Author: nasa   File: analysis.py    License: MIT License 6 votes vote down vote up
def _importShapefile(shapefile, dbname):
    """Import shapefile into database *dbname*."""
    db = dbio.connect(dbname)
    cur = db.cursor()
    tablename = ''.join(random.choice(string.ascii_lowercase) for _ in range(6))
    sql = "create table {0} (gid serial)".format(tablename)
    cur.execute(sql)
    cur.execute("select addgeometrycolumn('{0}', 'geom', 4326, 'POLYGON', 2)".format(tablename))
    ds = ogr.Open(shapefile)
    lyr = ds.GetLayer()
    for feat in lyr:
        geom = feat.GetGeometryRef()
        sql = "insert into {0} (geom) values (st_geometryfromtext('{1}', 4326))".format(tablename, geom.ExportToWkt())
        cur.execute(sql)
    ds.Destroy()
    db.commit()
    cur.close()
    db.close()
    return tablename 
Example #23
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 #24
Source Project: neon   Author: NervanaSystems   File: spacenet_utils.py    License: Apache License 2.0 6 votes vote down vote up
def load_as_uint8(filename):

    image = gdal.Open(filename)
    image_array = np.array(image.ReadAsArray())

    image_uint8 = np.zeros(image_array.shape, dtype=np.uint8)
    # rescale each band to be between 0, 255

    for k, band in enumerate(image_array):
        band_max = np.max(band)

        if band_max != 0:
            band = band.astype(np.float) / band_max * 255.0

        image_uint8[k, :, :] = band

    return image_uint8 
Example #25
Source Project: Carnets   Author: holzschu   File: test_shp.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_attributeexport(self):
        def testattributes(lyr, graph):
            feature = lyr.GetNextFeature()
            while feature:
                coords = []
                ref = feature.GetGeometryRef()
                last = ref.GetPointCount() - 1
                edge_nodes = (ref.GetPoint_2D(0), ref.GetPoint_2D(last))
                name = feature.GetFieldAsString('Name')
                assert_equal(graph.get_edge_data(*edge_nodes)['Name'], name)
                feature = lyr.GetNextFeature()

        tpath = os.path.join(tempfile.gettempdir(), 'shpdir')

        G = nx.read_shp(self.shppath)
        nx.write_shp(G, tpath)
        shpdir = ogr.Open(tpath)
        edges = shpdir.GetLayerByName("edges")
        testattributes(edges, G)

    # Test export of node attributes in nx.write_shp (#2778) 
Example #26
Source Project: Carnets   Author: holzschu   File: test_shp.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_wkt_export(self):
        G = nx.DiGraph()
        tpath = os.path.join(tempfile.gettempdir(), 'shpdir')
        points = (
            "POINT (0.9 0.9)",
            "POINT (4 2)"
        )
        line = (
            "LINESTRING (0.9 0.9,4 2)",
        )
        G.add_node(1, Wkt=points[0])
        G.add_node(2, Wkt=points[1])
        G.add_edge(1, 2, Wkt=line[0])
        try:
            nx.write_shp(G, tpath)
        except Exception as e:
            assert False, e
        shpdir = ogr.Open(tpath)
        self.checkgeom(shpdir.GetLayerByName("nodes"), points)
        self.checkgeom(shpdir.GetLayerByName("edges"), line) 
Example #27
Source Project: qgis-cartodb   Author: gkudos   File: CartoDBLayer.py    License: GNU General Public License v2.0 5 votes vote down vote up
def __init__(self, iface, tableName, user, apiKey, owner=None, sql=None, geoJSON=None,
                 filterByExtent=False, spatiaLite=None, readonly=False, multiuser=False, isSQL=False):
        self.iface = iface
        self.user = user
        self._apiKey = apiKey
        self.layerType = 'ogr'
        self.owner = owner
        self.isSQL = isSQL
        self.multiuser = multiuser
        # SQLite available?
        driverName = "SQLite"
        sqLiteDrv = ogr.GetDriverByName(driverName)
        self.database_path = QgisCartoDB.CartoDBPlugin.PLUGIN_DIR + '/db/database.sqlite'
        self.datasource = sqLiteDrv.Open(self.database_path, True)
        self.layerName = tableName
        self.cartoTable = tableName
        self.readonly = False
        self._deletedFeatures = []
        self.sql = sql

        self.forceReadOnly = False or readonly

        if sql is None:
            sql = 'SELECT * FROM ' + self._schema() + self.cartoTable
            if filterByExtent:
                extent = self.iface.mapCanvas().extent()
                sql = sql + " WHERE ST_Intersects(ST_GeometryFromText('{}', 4326), the_geom)".format(extent.asWktPolygon())
        else:
            self.forceReadOnly = True

        self._loadData(sql, geoJSON, spatiaLite) 
Example #28
Source Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_VectorTools.py    License: MIT License 5 votes vote down vote up
def readFile(filename):
	print("Hey buddy, Reading the file: "+filename)

	filehandle = gdal.Open(filename, GA_ReadOnly )
	if filehandle == None:
		raise Exception("Unable to read the data file")

	band1 = filehandle.GetRasterBand(1)
	geotransform = filehandle.GetGeoTransform()
	geoproj = filehandle.GetProjection()
	Z = band1.ReadAsArray()
	xsize = filehandle.RasterXSize
	ysize = filehandle.RasterYSize
	return xsize,ysize,geotransform,geoproj,Z 
Example #29
Source Project: LSDMappingTools   Author: LSDtopotools   File: LSD_GeologyTools.py    License: MIT License 5 votes vote down vote up
def readFile(filename):
    print("Hey buddy, Reading the file: "+filename)

    filehandle = gdal.Open(filename, GA_ReadOnly )
    if filehandle == None:
        raise Exception("Unable to read the data file")

    band1 = filehandle.GetRasterBand(1)
    geotransform = filehandle.GetGeoTransform()
    geoproj = filehandle.GetProjection()
    Z = band1.ReadAsArray()
    xsize = filehandle.RasterXSize
    ysize = filehandle.RasterYSize
    return xsize,ysize,geotransform,geoproj,Z 
Example #30
Source Project: pyeo   Author: clcr   File: raster_manipulation.py    License: GNU General Public License v3.0 5 votes vote down vote up
def reproject_image(in_raster, out_raster_path, new_projection,  driver = "GTiff",  memory = 2e3, do_post_resample=True):
    """
    Creates a new, reprojected image from in_raster using the gdal.ReprojectImage function.

    Parameters
    ----------
    in_raster
        Either a gdal.Raster object or a path to a raster
    out_raster_path
        The path to the new output raster.
    new_projection
        The new projection in .wkt
    driver
        The format of the output raster.
    memory
        The amount of memory to give to the reprojection.
    do_post_resample
        If set to false, do not resample the image back to the original projection.

    Notes
    -----
    The GDAL reprojection routine changes the size of the pixels by a very small amount; for example, a 10m pixel image
    can become a 10.002m pixel resolution image. To stop alignment issues,
    by deafult this function resamples the images back to their original resolution

    """
    if type(new_projection) is int:
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(new_projection)
        new_projection = proj.ExportToWkt()
    log = logging.getLogger(__name__)
    log.info("Reprojecting {} to {}".format(in_raster, new_projection))
    if type(in_raster) is str:
        in_raster = gdal.Open(in_raster)
    res = in_raster.GetGeoTransform()[1]
    gdal.Warp(out_raster_path, in_raster, dstSRS=new_projection, warpMemoryLimit=memory, format=driver)
    # After warping, image has irregular gt; resample back to previous pixel size
    # TODO: Make this an option
    if do_post_resample:
        resample_image_in_place(out_raster_path, res)
    return out_raster_path