Python osgeo.ogr.GetDriverByName() Examples

The following are code examples for showing how to use osgeo.ogr.GetDriverByName(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: lidar   Author: giswqs   File: slicing.py    MIT License 7 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)
    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# convert images in a selected folder to shapefiles 
Example 2
Project: core   Author: lifemapper   File: vector_validator.py    GNU General Public License v3.0 6 votes vote down vote up
def validate_vector_file(vector_filename):
    """Validates a vector file by seeing if it can be loaded by OGR

    Args:
        vector_filename : The file path of a vector file to validate
    """
    msg = 'Valid'
    valid = False
    if os.path.exists(vector_filename):
        try:
            drv = ogr.GetDriverByName(LMFormat.SHAPE.driver)
            ds = drv.Open(vector_filename)
            if ds is None:
                msg = 'Could not open {}'.format(vector_filename)
            else:
                lyr = ds.GetLayer()
                feature_count = lyr.GetFeatureCount()
                valid = True
        except Exception, e:
            msg = str(e) 
Example 3
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def getShapefileRowHeaders(cls, shapefileFilename):
        """
        @summary: Get a (sorted by feature id) list of row headers for a shapefile
        @todo: Move this to a common Vector class in LmCommon or LmBackend and 
                 use with LmCompute.  This is a rough copy of what is now used for 
                 rasterIntersect.
        """
        ogr.RegisterAll()
        drv = ogr.GetDriverByName(LMFormat.getDefaultOGR().driver)
        ds = drv.Open(shapefileFilename)
        lyr = ds.GetLayer(0)
        
        rowHeaders = []
        
        for j in range(lyr.GetFeatureCount()):
            curFeat = lyr.GetFeature(j)
            siteIdx = curFeat.GetFID()
            x, y = curFeat.geometry().Centroid().GetPoint_2D()
            rowHeaders.append((siteIdx, x, y))
            
        return sorted(rowHeaders)

# ............................................... 
Example 4
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def copyData(self, sourceDataLocation, targetDataLocation=None, 
                     format=LMFormat.getDefaultOGR().driver):
        """
        Copy sourceDataLocation dataset to targetDataLocation or this layer's 
        dlocation.
        """
        if sourceDataLocation is not None and os.path.exists(sourceDataLocation):
            if targetDataLocation is None:
                if self._dlocation is not None:
                    targetDataLocation = self._dlocation
                else:
                    raise LMError('Target location is None')
        else:
            raise LMError('Source location %s is invalid' % str(sourceDataLocation))

        ogr.RegisterAll()
        drv = ogr.GetDriverByName(format)
        try:
            ds = drv.Open(sourceDataLocation)
        except Exception, e:
            raise LMError(['Invalid datasource' % sourceDataLocation, str(e)]) 
Example 5
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def testVector(dlocation, driver=LMFormat.getDefaultOGR().driver):
        goodData = False
        featCount = 0
        if dlocation is not None and os.path.exists(dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(driver)
            try:
                ds = drv.Open(dlocation)
            except Exception, e:
                goodData = False
            else:
                try:
                    slyr = ds.GetLayer(0)
                except Exception, e:
                    goodData = False
                else: 
Example 6
Project: core   Author: lifemapper   File: assemble_pam_from_solr.py    GNU General Public License v3.0 6 votes vote down vote up
def _get_row_headers_from_shapegrid(shapegrid_filename):
    """Get the row headers from a shapegrid

    Args:
        shapegrid_filename (:obj: `str`): The file location of the shapegrid to
            get headers
    """
    ogr.RegisterAll()
    drv = ogr.GetDriverByName(LMFormat.getDefaultOGR().driver)
    ds = drv.Open(shapegrid_filename)
    lyr = ds.GetLayer(0)
    
    row_headers = []
    
    for j in range(lyr.GetFeatureCount()):
        curFeat = lyr.GetFeature(j)
        siteIdx = curFeat.GetFID()
        x, y = curFeat.geometry().Centroid().GetPoint_2D()
        row_headers.append((siteIdx, x, y))
    return sorted(row_headers)

# ............................................................................. 
Example 7
Project: core   Author: lifemapper   File: subset.py    GNU General Public License v3.0 6 votes vote down vote up
def getRowHeaders(shapefileFilename):
    """
    @summary: Get a (sorted by feature id) list of row headers for a shapefile
    @todo: Move this to LmCommon or LmBackend and use with LmCompute.  This is
                 a rough copy of what is now used for rasterIntersect
    """
    ogr.RegisterAll()
    drv = ogr.GetDriverByName(LMFormat.getDefaultOGR().driver)
    ds = drv.Open(shapefileFilename)
    lyr = ds.GetLayer(0)
    
    rowHeaders = []
    
    for j in range(lyr.GetFeatureCount()):
        curFeat = lyr.GetFeature(j)
        siteIdx = curFeat.GetFID()
        x, y = curFeat.geometry().Centroid().GetPoint_2D()
        rowHeaders.append((siteIdx, x, y))
        
    return sorted(rowHeaders) 
Example 8
Project: core   Author: lifemapper   File: createshape.py    GNU General Public License v3.0 6 votes vote down vote up
def testShapefile(dlocation):
        """
        @todo: This should go into a LmCommon base layer class
        """
        goodData = True
        featCount = 0
        if dlocation is not None and os.path.exists(dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(LMFormat.getDefaultOGR().driver)
            try:
                ds = drv.Open(dlocation)
            except Exception, e:
                goodData = False
            else:
                try:
                    slyr = ds.GetLayer(0)
                except Exception, e:
                    goodData = False
                else: 
Example 9
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_VectorTools.py    MIT License 6 votes vote down vote up
def writeFile(filename,geotransform,geoprojection,data):
	(x,y) = data.shape
	format = "GTiff"
	noDataValue = -9999
	driver = gdal.GetDriverByName(format)
	# you can change the dataformat but be sure to be able to store negative values including -9999
	dst_datatype = gdal.GDT_Float32

	#print(data)

	dst_ds = driver.Create(filename,y,x,1,dst_datatype)
	dst_ds.GetRasterBand(1).WriteArray(data)
	dst_ds.GetRasterBand(1).SetNoDataValue( noDataValue )
	dst_ds.SetGeoTransform(geotransform)
	dst_ds.SetProjection(geoprojection)
	return 1 
Example 10
Project: LSDMappingTools   Author: LSDtopotools   File: LSD_GeologyTools.py    MIT License 6 votes vote down vote up
def writeFile(filename,geotransform,geoprojection,data):
    (x,y) = data.shape
    format = "GTiff"
    noDataValue = -9999
    driver = gdal.GetDriverByName(format)
    # you can change the dataformat but be sure to be able to store negative values including -9999
    dst_datatype = gdal.GDT_Float32

    #print(data)

    dst_ds = driver.Create(filename,y,x,1,dst_datatype)
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds.GetRasterBand(1).SetNoDataValue( noDataValue )
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.SetProjection(geoprojection)
    return 1 
Example 11
Project: lidar   Author: giswqs   File: slicing.py    MIT License 6 votes vote down vote up
def writeRaster(arr, out_path, template):
    no_data = 0
    # First of all, gather some information from the template file
    data = gdal.Open(template)
    [cols, rows] = arr.shape
    trans = data.GetGeoTransform()
    proj = data.GetProjection()
    # nodatav = 0 #data.GetNoDataValue()
    # Create the file, using the information from the template file
    outdriver = gdal.GetDriverByName("GTiff")
    # http://www.gdal.org/gdal_8h.html
    # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,
    outdata   = outdriver.Create(str(out_path), rows, cols, 1, gdal.GDT_UInt32)
    # Write the array to the file, which is the original array in this example
    outdata.GetRasterBand(1).WriteArray(arr)
    # Set a no data value if required
    outdata.GetRasterBand(1).SetNoDataValue(no_data)
    # Georeference the image
    outdata.SetGeoTransform(trans)
    # Write projection information
    outdata.SetProjection(proj)
    return arr


# raster to vector 
Example 12
Project: unmixing   Author: arthur-e   File: lsma.py    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 13
Project: Planet-GEE-Pipeline-GUI   Author: samapriya   File: kml_aoi.py    Apache License 2.0 6 votes vote down vote up
def parsed(args):
        kml_file=args.geo
        def kml2geojson(kml_file):
            drv = ogr.GetDriverByName('KML')
            kml_ds = drv.Open(kml_file)
            for kml_lyr in kml_ds:
                for feat in kml_lyr:
                    outfile=feat.ExportToJson()
                    geom2=str(outfile).replace(", 0.0",'')
                    with open(args.loc+'./kmlout.geojson','w') as csvfile:
                        writer=csv.writer(csvfile)
                        writer.writerow([geom2])
        kml2geojson(args.geo)
        raw= open(args.loc+'./kmlout.geojson')
        for line in raw:
            fields=line.strip().split(":")[3]
            f2=fields.strip().split("}")[0]
            filenames = p1+f2+p2+str(args.start)+p3+str(args.end)+p4+p5+str(args.cloud)+p6
        with open(args.loc+'./aoi.json', 'w') as outfile:
            outfile.write(filenames)
            outfile.close() 
Example 14
Project: Bandit   Author: paknorton   File: prms_geo.py    MIT License 6 votes vote down vote up
def __init__(self, filename):
        """Create the Geo object.

        :param str filename: name of file
        """

        err = GdalErrorHandler()
        handler = err.handler

        gdal.PushErrorHandler(handler)
        gdal.UseExceptions()

        self.__filename = filename
        self.__layers = None
        self.__selected_layer = None

        # use OGR specific exceptions
        ogr.UseExceptions()

        # Load the file - this assumes a file geodatabase
        driver = ogr.GetDriverByName('OpenFileGDB')
        self.__gdb = driver.Open(self.__filename) 
Example 15
Project: eetc   Author: george-azzari   File: gdaltools.py    MIT License 6 votes vote down vote up
def get_polygon_corners(self, affine=True, **polyargs):
        """ Returns the extent of the polygon in georeferenced or affine coordinates.
            Georeferenced coordinates are directly read from the polygon file.
            Affine coordinates are calculated using the transformation set with the current instance of the class."""
        d = ogr.GetDriverByName(polyargs['driver'])
        poly = d.Open(polyargs['fpath'])
        layer = poly.GetLayerByName(polyargs['layer_name'])
        ext = layer.GetExtent()
        ulx, lrx, lry, uly = ext[0], ext[1], ext[2], ext[3]
        if not affine:
            return ulx, lrx, lry, uly
        if affine:
            lin_ul = self.get_affinecoord(ulx, uly)
            lin_lr = self.get_affinecoord(lrx, lry)
            top, bottom, left, right = lin_ul[1], lin_lr[1], lin_ul[0], lin_lr[0]
            return left, bottom, right, top  # this is the order for warping 
Example 16
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions.py    MIT License 6 votes vote down vote up
def toRaster(array, path, geometry, srs, navalue=-9999):
    """
    Writes a single array to a raster with coordinate system and geometric
    information.

    path = target path
    srs = spatial reference system
    """
    xpixels = array.shape[1]
    ypixels = array.shape[0]
    path = path.encode('utf-8')
    image = gdal.GetDriverByName("GTiff").Create(path, xpixels, ypixels,
                                1, gdal.GDT_Float32)
    image.SetGeoTransform(geometry)
    image.SetProjection(srs)
    image.GetRasterBand(1).WriteArray(array)
    image.GetRasterBand(1).SetNoDataValue(navalue) 
Example 17
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions.py    MIT License 6 votes vote down vote up
def toRasters(arraylist, path, geometry, srs):
    """
    Writes a list of 2d arrays, or a 3d array, to a series of tif files.

    Arraylist format = [[name,array],[name,array],....]
    path = target path
    geometry = gdal geometry object
    srs = spatial reference system object
    """
    if path[-2:] == "\\":
        path = path
    else:
        path = path + "\\"
    sample = arraylist[0][1]
    ypixels = sample.shape[0]
    xpixels = sample.shape[1]
    for ray in  tqdm(arraylist):
        image = gdal.GetDriverByName("GTiff").Create(os.path.join(path,
                                                              ray[0] + ".tif"),
                                    xpixels, ypixels, 1, gdal.GDT_Float32)
        image.SetGeoTransform(geometry)
        image.SetProjection(srs)
        image.GetRasterBand(1).WriteArray(ray[1]) 
Example 18
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions_v3.py    MIT License 6 votes vote down vote up
def toRaster(array, path, geometry, srs, navalue=-9999):
    """
    Writes a single array to a raster with coordinate system and geometric
    information.

    path = target path
    srs = spatial reference system
    """
    xpixels = array.shape[1]    
    ypixels = array.shape[0]
    path = path.encode('utf-8')
    image = gdal.GetDriverByName("GTiff").Create(path, xpixels, ypixels,
                                1, gdal.GDT_Float32)
    image.SetGeoTransform(geometry)
    image.SetProjection(srs)
    image.GetRasterBand(1).WriteArray(array)
    image.GetRasterBand(1).SetNoDataValue(navalue) 
Example 19
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions_v3.py    MIT License 6 votes vote down vote up
def toRasters(arraylist, path, geometry, srs):
    """
    Writes a list of 2d arrays, or a 3d array, to a series of tif files.

    Arraylist format = [[name,array],[name,array],....]
    path = target path
    geometry = gdal geometry object
    srs = spatial reference system object
    """
    if path[-2:] == "\\":
        path = path
    else:
        path = path + "\\"
    sample = arraylist[0][1]
    ypixels = sample.shape[0]
    xpixels = sample.shape[1]
    for ray in  tqdm(arraylist):
        image = gdal.GetDriverByName("GTiff").Create(os.path.join(path,
                                                              ray[0] + ".tif"),
                                    xpixels, ypixels, 1, gdal.GDT_Float32)
        image.SetGeoTransform(geometry)
        image.SetProjection(srs)
        image.GetRasterBand(1).WriteArray(ray[1]) 
Example 20
Project: PyGeoC   Author: lreis2415   File: vector.py    MIT License 6 votes vote down vote up
def write_line_shp(line_list, out_shp):
        """Export ESRI Shapefile -- Line feature"""
        print('Write line shapefile: %s' % out_shp)
        driver = ogr_GetDriverByName(str('ESRI Shapefile'))
        if driver is None:
            print('ESRI Shapefile driver not available.')
            sys.exit(1)
        if os.path.exists(out_shp):
            driver.DeleteDataSource(out_shp)
        ds = driver.CreateDataSource(out_shp.rpartition(os.sep)[0])
        if ds is None:
            print('ERROR Output: Creation of output file failed.')
            sys.exit(1)
        lyr = ds.CreateLayer(str(out_shp.rpartition(os.sep)[2].split('.')[0]), None, wkbLineString)
        for l in line_list:
            line = ogr_Geometry(wkbLineString)
            for i in l:
                line.AddPoint(i[0], i[1])
            templine = ogr_CreateGeometryFromJson(line.ExportToJson())
            feature = ogr_Feature(lyr.GetLayerDefn())
            feature.SetGeometry(templine)
            lyr.CreateFeature(feature)
            feature.Destroy()
        ds.Destroy() 
Example 21
Project: coast-def   Author: zdb999   File: utils.py    GNU General Public License v3.0 6 votes vote down vote up
def reproject(source_path, match_path, out_path ):
  # Source
  src_filename = real_path + source_path
  src = gdal.Open(src_filename, gdalconst.GA_ReadOnly)
  src_proj = src.GetProjection()
  src_geotrans = src.GetGeoTransform()

  # We want a section of source that matches this:
  match_filename = real_path + match_path
  match_ds = gdal.Open(match_filename, gdalconst.GA_ReadOnly)
  match_proj = match_ds.GetProjection()
  match_geotrans = match_ds.GetGeoTransform()
  wide = match_ds.RasterXSize
  high = match_ds.RasterYSize

  # Output / destination
  dst_filename = real_path + out_path
  dst = gdal.GetDriverByName('GTiff').Create(dst_filename, wide, high, 1, gdalconst.GDT_Float32)
  dst.SetGeoTransform(match_geotrans)
  dst.SetProjection( match_proj)

  # Do the work
  gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)

  del dst # Flush 
Example 22
Project: core   Author: lifemapper   File: parameter_sweep.py    GNU General Public License v3.0 5 votes vote down vote up
def _get_model_points(self, occ_shp_filename):
        """Get minimal point csv to be used for modeling.

        Args:
            occ_shp_filename : The file location of a point shapefile.
        
        Note:
            * Removes duplicate point locations
        """
        points = set([])
        drv = ogr.GetDriverByName(LMFormat.SHAPE.driver)
        ds = drv.Open(occ_shp_filename, 0)
        lyr = ds.GetLayer()

        for feature in lyr:
            geom = feature.GetGeometryRef()
            pt_geom = geom.GetPoint()
            points.add((pt_geom[0], pt_geom[1]))
            
        # Add identifiers and create list
        i = 0
        ret_points = []
        for x, y in points:
            ret_points.append((i, x, y))
            i += 1

        return ret_points

    # ........................................ 
Example 23
Project: core   Author: lifemapper   File: create_mask.py    GNU General Public License v3.0 5 votes vote down vote up
def write_tiff(out_filename, bbox, cell_size, data, epsg,
               nodata=DEFAULT_NODATA):
    """
    @summary: Write a GeoTiff raster layer from a numpy array
    @param out_filename: The file location to write the raster to
    @param bbox: (minx, miny, maxx, maxy) tuple of raster coordinates
    @param cell_size: The cell size, in map units, of each cell in the raster
    @param data: The data array to use for raster content
    @param epsg: The epsg code of the map projection to use for this raster
    @param nodata: A value to use for NODATA
    @todo: Probably should establish a common layer package in LmCommon
    """
    minx, miny, maxx, maxy = bbox
    (num_rows, num_cols) = data.shape
    
    readyFilename(out_filename, overwrite=True)
    
    drv = gdal.GetDriverByName(LMFormat.GTIFF.driver)
    ds = drv.Create(out_filename, num_cols, num_rows, 1, gdalconst.GDT_Byte)
    ds.SetGeoTransform([minx, cell_size, 0, maxy, 0, -cell_size])

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    ds.SetProjection(srs.ExportToWkt())
    
    outBand = ds.GetRasterBand(1)
    outBand.WriteArray(data)
    outBand.FlushCache()
    outBand.SetNoDataValue(nodata)

# .............................................................................
# Testing code 
Example 24
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def getFormatLongName(self):
        name = ''
        drv = gdal.GetDriverByName(self._dataFormat)
        if drv is not None:
            name = drv.GetMetadataItem('DMD_LONGNAME')
        return name
            
# ............................................... 
Example 25
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def _copyGDALData(self, bandnum, infname, outfname, format='GTiff', kwargs={}):
        """
        @summary: Copy the dataset into a new file.  
        @param bandnum: The band number to read.
        @param outfname: Filename to write this dataset to.
        @param format: GDAL-writeable raster format to use for new dataset. 
                            http://www.gdal.org/formats_list.html
        @param doFlip: True if data begins at the southern edge of the region
        @param doShift: True if the leftmost edge of the data should be shifted 
                 to the center (and right half shifted around to the beginning) 
        @param nodata: Value used to indicate nodata in the new file.
        @param srs: Spatial reference system to use for the data. This is only 
                        necessary if the dataset does not have an SRS present.  This
                        will NOT project the dataset into a different projection.
        """
        options = []
        if format == 'AAIGrid':
            options = ['FORCE_CELLSIZE=True']
            #kwargs['FORCE_CELLSIZE'] = True
            #kwargs['DECIMAL_PRECISION'] = 4
        driver = gdal.GetDriverByName(format)
        metadata = driver.GetMetadata()
        if not (metadata.has_key(gdal.DCAP_CREATECOPY) 
                     and metadata[gdal.DCAP_CREATECOPY] == 'YES'):
            raise LMError(currargs='Driver %s does not support CreateCopy() method.' 
                              % format)
        inds = gdal.Open( infname )
        try:
            outds = driver.CreateCopy(outfname, inds, 0, options)
        except Exception, e:
            raise LMError(currargs='Creation failed for %s from band %d of %s (%s)'
                                          % (outfname, bandnum, infname, str(e))) 
Example 26
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def verifyField(self, dlocation, ogrFormat, attrname):
        """
        @summary: Read OGR-accessible data and save the features and 
                     featureAttributes on the Vector object
        @param dlocation: Location of the data
        @param ogrFormat: OGR-supported data format code, available at
                             http://www.gdal.org/ogr/ogr_formats.html
        @return: boolean for success/failure 
        @raise LMError: on failure to read data.
        @note: populateStats calls this
        """
        if attrname is None:
            fieldOk = True
        else:
            fieldOk = False
            if dlocation is not None and os.path.exists(dlocation):
                ogr.RegisterAll()
                drv = ogr.GetDriverByName(ogrFormat)
                try:
                    ds = drv.Open(dlocation)
                except Exception, e:
                    raise LMError(['Invalid datasource' % dlocation, str(e)])
                
                lyrDef = ds.GetLayer(0).GetLayerDefn()
                # Make sure given field exists and is the correct type
                for i in range(lyrDef.GetFieldCount()):
                    fld = lyrDef.GetFieldDefn(i)
                    fldname = fld.GetNameRef()
                    if attrname == fldname:
                        fldtype = fld.GetType()
                        if fldtype in (ogr.OFTInteger, ogr.OFTReal, ogr.OFTBinary):
                            fieldOk = True
                        break 
Example 27
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def readWithOGR(self, dlocation, ogrFormat, featureLimit=None, doReadData=False):
        """
        @summary: Read OGR-accessible data and set the features and 
                     featureAttributes on the Vector object
        @param dlocation: Location of the data
        @param ogrFormat: OGR-supported data format code, available at
                             http://www.gdal.org/ogr/ogr_formats.html
        @return: boolean for success/failure 
        @raise LMError: on failure to read data.
        @note: populateStats calls this
        @todo: remove featureLimit, read subsetDLocation if there is a limit 
        """
        thisBBox = localIdIdx = geomIdx = feats = featAttrs = None
        if dlocation is not None and os.path.exists(dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(ogrFormat)
            try:
                ds = drv.Open(dlocation)
            except Exception, e:
                raise LMError(['Invalid datasource' % dlocation, str(e)])
                          
            self.clearFeatures() 
            try:
                slyr = ds.GetLayer(0)
            except Exception, e:
                raise LMError(currargs='#### Failed to GetLayer from %s' % dlocation,
                                  prevargs=e.args, doTrace=True)

            # Get bounding box 
Example 28
Project: core   Author: lifemapper   File: create_mask_tiff.py    GNU General Public License v3.0 5 votes vote down vote up
def getXYsForShapefile(filename):
   def getXY(wkt):
      startidx = wkt.find('(')
      if wkt[:startidx].strip().lower() == 'point':
         tmp = wkt[startidx+1:]
         endidx = tmp.find(')')
         tmp = tmp[:endidx]
         vals = tmp.split()
         if len(vals)  == 2:
            try:
               x = float(vals[0])
               y = float(vals[1])
               return x, y
            except:
               return None
         else:
            return None
      else:
         return None

   drv = ogr.GetDriverByName('ESRI Shapefile')
   ds = drv.Open(filename)
   lyr = ds.GetLayer()
   
   pts = []
   for feat in lyr:
      geom = feat.GetGeometryRef()
      r = getXY(geom.Centroid().ExportToWkt())
      if r is not None:
         pts.append(r)
   lyr = ds = None
   
   return pts 
Example 29
Project: HAPI   Author: MAfarrag   File: GISpy.py    MIT License 5 votes vote down vote up
def SaveRaster(raster,path):
    """
    ===================================================================
      SaveRaster(raster,path)
    ===================================================================
    this function saves a raster to a path
    
    inputs:
    ----------
        1- raster:
            [gdal object]
        2- path:
            [string] a path includng the name of the raster and extention like 
            path="data/cropped.tif"
    
    Outputs:
    ----------
        the function does not return and data but only save the raster to the hard drive
    
    EX:
    ----------
        SaveRaster(raster,output_path)
    """
    #### input data validation
    # data type
    assert type(raster)==gdal.Dataset, "src should be read using gdal (gdal dataset please read it using gdal library) "
    assert type(path)== str, "Raster_path input should be string type"
    # input values
    ext=path[-4:]
    assert ext == ".tif", "please add the extension at the end of the path input"
    
    driver = gdal.GetDriverByName ( "GTiff" )
    dst_ds = driver.CreateCopy( path, raster, 0 )
    dst_ds = None # Flush the dataset to disk 
Example 30
Project: HAPI   Author: MAfarrag   File: GISpy.py    MIT License 5 votes vote down vote up
def SaveRaster(raster,path):
    """
    ===================================================================
      SaveRaster(raster,path)
    ===================================================================
    this function saves a raster to a path
    
    inputs:
    ----------
        1- raster:
            [gdal object]
        2- path:
            [string] a path includng the name of the raster and extention like 
            path="data/cropped.tif"
    
    Outputs:
    ----------
        the function does not return and data but only save the raster to the hard drive
    
    EX:
    ----------
        SaveRaster(raster,output_path)
    """
    #### input data validation
    # data type
    assert type(raster)==gdal.Dataset, "src should be read using gdal (gdal dataset please read it using gdal library) "
    assert type(path)== str, "Raster_path input should be string type"
    # input values
    ext=path[-4:]
    assert ext == ".tif", "please add the extension at the end of the path input"
    
    driver = gdal.GetDriverByName ( "GTiff" )
    dst_ds = driver.CreateCopy( path, raster, 0 )
    dst_ds = None # Flush the dataset to disk 
Example 31
Project: wikiwhere   Author: mkrnr   File: countries.py    MIT License 5 votes vote down vote up
def __init__(self, country_file):
        driver = ogr.GetDriverByName('ESRI Shapefile')
        self.countryFile = driver.Open(country_file)
        self.layer = self.countryFile.GetLayer() 
Example 32
Project: Data-Tools-for-ArcGIS   Author: datastory   File: Service Geometry Downloader.py    MIT License 5 votes vote down vote up
def join_files(folder):
    """
    Merge all files in given folder into one datasource

    :param folder: name of temporary folder to be joined
    :return: resulting joined files file name
    """

    first_ds = _get_first_source_dataset(folder)
    first_layer = first_ds.GetLayer(0)

    drv = ogr.GetDriverByName('GeoJSON')
    tmpfile = tempfile.mktemp(suffix=".json")
    out_dsn = drv.CreateDataSource(tmpfile)
    out_layer = out_dsn.CopyLayer(first_layer, new_name=first_layer.GetName())

    for source in os.listdir(folder)[1:]:
        logging.info('Joining file %s to %s' % (source, tmpfile)) 
        dsn = ogr.Open(os.path.join(folder, source))
        layer = dsn.GetLayer()
        nfeatures = layer.GetFeatureCount()

        for i in range(nfeatures):
            feature = layer.GetNextFeature()
            out_layer.CreateFeature(feature.Clone())

    out_dsn.Destroy()

    return tmpfile 
Example 33
Project: qgis-geogiglight-plugin   Author: planetfederal   File: testgpkg.py    GNU General Public License v2.0 5 votes vote down vote up
def _getOgrTestLayer(self):
        dest = self._copyTestLayer()
        driver = ogr.GetDriverByName('GPKG')
        dataSource = driver.Open(dest, 1)
        return dataSource 
Example 34
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    MIT License 5 votes vote down vote up
def array2raster(rasterfn,newRasterfn,array,driver_name = "ENVI", noDataValue = -9999):
    """Takes an array and writes to a GDAL compatible raster. It needs another raster to map the dimensions.

    Args:
        FileName (str): The filename (with path and extension) of a raster that has the same dimensions as the raster to be written.
        newRasterfn (str): The filename (with path and extension) of the new raster.
        array (np.array): The array to be written
        driver_name (str): The type of raster to write. Default is ENVI since that is the LSDTOpoTools format
        noDataValue (float): The no data value

    Return:
        np.array: A numpy array with the data from the raster.

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName(driver_name)
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outRaster.GetRasterBand(1).SetNoDataValue( noDataValue )
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
#============================================================================== 
Example 35
Project: lidar   Author: giswqs   File: filling.py    MIT License 5 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)

    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    # raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# extract sinks from dem 
Example 36
Project: coded   Author: bullocke   File: classify.py    MIT License 5 votes vote down vote up
def create_mask_from_vector(vector_data_path, cols, rows, geo_transform, 
			    projection, target_value=1,
                            output_fname='', dataset_format='MEM'):

    """
    Rasterize the given vector (wrapper for gdal.RasterizeLayer). 
    Return a gdal.Dataset.
    :param vector_data_path: Path to a shapefile
    :param cols: Number of columns of the result
    :param rows: Number of rows of the result
    :param geo_transform: Returned value of gdal.Dataset.GetGeoTransform 
	(coefficients for transforming between pixel/line (P,L) raster space,
	 and projection coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by 
	gdal.Dataset.GetProjectionRef)
    :param target_value: Pixel value for the pixels. Must be a valid 
	gdal.GDT_UInt16 value.
    :param output_fname: If the dataset_format is GeoTIFF, this is the output 
	file name
    :param dataset_format: The gdal.Dataset driver name. [default: MEM]
    """

    driver = ogr.GetDriverByName('ESRI Shapefile')
    data_source = driver.Open(vector_data_path, 0)
    if data_source is None:
        report_and_exit("File read failed: %s", vector_data_path)
    layer = data_source.GetLayer(0)
    driver = gdal.GetDriverByName(dataset_format)
    target_ds = driver.Create(output_fname, cols, rows, 1, gdal.GDT_UInt16)
    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(projection)
    gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[target_value])
    return target_ds 
Example 37
Project: coded   Author: bullocke   File: classify.py    MIT License 5 votes vote down vote up
def write_geotiff(fname, data, geo_transform, projection, data_type=gdal.GDT_Byte):

    """
    Create a GeoTIFF file with the given data.
    :param fname: Path to a directory with shapefiles
    :param data: Number of rows of the result
    :param geo_transform: Returned value of 
	gdal.Dataset.GetGeoTransform (coefficients for transforming between 
	pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by 
	gdal.Dataset.GetProjectionRef)
    """

    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    dataset = driver.Create(fname, cols, rows, 1, data_type)
    dataset.SetGeoTransform(geo_transform)
    dataset.SetProjection(projection)
    band = dataset.GetRasterBand(1)
    band.WriteArray(data)

    metadata = {
        'TIFFTAG_COPYRIGHT': 'CC BY 4.0',
        'TIFFTAG_DOCUMENTNAME': 'classification',
        'TIFFTAG_IMAGEDESCRIPTION': 'Supervised classification.',
        'TIFFTAG_SOFTWARE': 'Python, GDAL, scikit-learn'
    }
    dataset.SetMetadata(metadata)

    dataset = None  # Close the file
    return 
Example 38
Project: coded   Author: bullocke   File: train_deg.py    MIT License 5 votes vote down vote up
def make_class_dict(path):
    # Set up dict to save Xs and Ys
    driver = ogr.GetDriverByName('ESRI Shapefile')
#    data_source = gdal.OpenEx(vector_data_path, gdal.OF_VECTOR)
    data_source = driver.Open(path, 0)
    if data_source is None:
        report_and_exit("File read failed: %s", path)

    layer = data_source.GetLayer(0)
    class_labels = []
    data = []

    for feature in layer:
        try:
            var1 = float(feature.GetField('NFDI_mag')) 
#   	    var2 = float(feature.GetField('NFDI_rmse'))
	    var3 = float(feature.GetField('NFDI_sin')) 
	    var4 = float(feature.GetField('NFDI_cos')) 
	    var5 = float(feature.GetField('Gv_mag')) 
	    var6 = float(feature.GetField('Shade_mag'))
	    var7 = float(feature.GetField('NPV_mag')) 
	    var8 = float(feature.GetField('Soil_mag')) 
	    label = feature.GetField('class')
        except:
            continue 
        class_labels.append(label)
        #data.append([var1, var2, var3, var4, var5, var6, var7, var8])
        data.append([var1, var3, var4, var5, var6, var7, var8])

    return class_labels, data 
Example 39
Project: coded   Author: bullocke   File: train.py    MIT License 5 votes vote down vote up
def make_class_dict(path):
    # Set up dict to save Xs and Ys
    driver = ogr.GetDriverByName('ESRI Shapefile')
    data_source = driver.Open(path, 0)
    if data_source is None:
        report_and_exit("File read failed: %s", path)

    layer = data_source.GetLayer(0)
    class_labels = []
    data = []

    for feature in layer:
        try:
            var1 = float(feature.GetField('NFDI_mag')) 
   	    var2 = float(feature.GetField('NFDI_rmse'))
	    var3 = float(feature.GetField('NFDI_sin')) 
	    var4 = float(feature.GetField('NFDI_cos')) 
	    var5 = float(feature.GetField('Gv_mag')) 
	    var6 = float(feature.GetField('Shade_mag'))
	    var7 = float(feature.GetField('NPV_mag')) 
	    var8 = float(feature.GetField('Soil_mag')) 
	    label = feature.GetField('class')
        except:
            continue 
        class_labels.append(label)
        data.append([var1, var2, var3, var4, var5, var6, var7, var8])
#        data.append([var1, var3, var4, var5, var6, var7, var8])

    return class_labels, data 
Example 40
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    GNU General Public License v3.0 5 votes vote down vote up
def write_geometry(geometry, out_path, srs_id=4326):
    """
    Saves the geometry in an ogr.Geometry object to a shapefile.

    Parameters
    ----------
    geometry
        An ogr.Geometry object
    out_path
        The location to save the output shapefile
    srs_id
        The projection of the output shapefile. Can be an EPSG number or a WKT string.

    Notes
    -----
    The shapefile consists of one layer named 'geometry'.


    """
    # TODO: Fix this needing an extra filepath on the end
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    if type(srs_id) is int:
        srs.ImportFromEPSG(srs_id)
    if type(srs_id) is str:
        srs.ImportFromWkt(srs_id)
    layer = data_source.CreateLayer(
        "geometry",
        srs,
        geom_type=geometry.GetGeometryType())
    feature_def = layer.GetLayerDefn()
    feature = ogr.Feature(feature_def)
    feature.SetGeometry(geometry)
    layer.CreateFeature(feature)
    data_source.FlushCache()
    data_source = None 
Example 41
Project: swat-tools   Author: DigitalAgricultureDiscovery   File: geotools.py    GNU General Public License v3.0 5 votes vote down vote up
def get_shapefile_extent(shapefile_filepath):
    """
    Opens up shapefile with GDAL/OGR and returns the
    bounding box (extent) for the shapefile.

    Parameters
    ----------
    shapefile_filepath: string
        location of the uploaded shapefile

    Returns
    -------
    extent: list
        list containing the extent (bounding box) for the shapefile

    """
    # load ESRI shapefile driver
    driver = ogr.GetDriverByName('ESRI Shapefile')

    # open the shapefile
    datasource = driver.Open(shapefile_filepath, 0)

    # get layer from shapefile
    layer = datasource.GetLayer()

    # get shapefile's extent
    extent = layer.GetExtent()

    return extent 
Example 42
Project: HackZurich   Author: RefugeeMatchmaking   File: test_shp.py    MIT License 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 43
Project: surfclass   Author: Kortforsyningen   File: vectorize.py    MIT License 5 votes vote down vote up
def open_or_create_destination_datasource(dst_ds_name, dst_format=None, dsco=None):
    """Open an OGR DataSource for update if it exists. Otherwise create it.

    Args:
        dst_ds_name (str): DataSource string.
        dst_format (str, optional): DataSource format as specified by OGR format string. Defaults to None.
        dsco (list of str, optional): List of OGR DataSource creation options. Defaults to None.

    Raises:
        Exception: If output datasource exists, but cannot be opened in update mode.
        Exception: If `dst_format` does not identify an OGR driver.
        Exception: Datasource cannot be created.

    Returns:
        osgeo.ogr.DataSource: OGR datasource.

    """
    logger.debug("Try to open '%s' for update", dst_ds_name)
    dst_ds = ogr.Open(dst_ds_name, update=1)
    if dst_ds is None:
        dst_ds = ogr.Open(dst_ds_name)
        if dst_ds is not None:
            raise Exception(
                'Output datasource "%s" exists, but cannot be opened in update mode'
                % dst_ds
            )
        drv = ogr.GetDriverByName(dst_format)
        if drv is None:
            raise Exception("Cannot find driver %s" % dst_format)
        logger.debug("Try to create '%s'", dst_ds_name)
        dst_ds = drv.CreateDataSource(dst_ds_name, options=dsco or [])
        if dst_ds is None:
            raise Exception("Cannot create datasource '%s'" % dst_ds_name)
    return dst_ds 
Example 44
Project: surfclass   Author: Kortforsyningen   File: test_vectors.py    MIT License 5 votes vote down vote up
def test_classcounter(classraster_filepath, polygons_filepath):
    vecreader = FeatureReader(polygons_filepath)
    rasreader = MaskedRasterReader(classraster_filepath)
    # Create output ds
    mem_drv = ogr.GetDriverByName("Memory")
    out_ds = mem_drv.CreateDataSource("out")
    out_lyr = open_or_create_similar_layer(vecreader.lyr, out_ds)

    in_features = list(vecreader)
    vecreader.reset_reading()

    classes = range(6)
    calc = ClassCounter(vecreader, rasreader, out_lyr, classes)
    calc.process()

    out_reader = FeatureReader(out_ds, out_lyr)
    out_reader.reset_reading()

    out_features = list(out_reader)
    assert len(out_features) == len(in_features)
    out_fields = set(out_features[0].keys())
    expected_classes = ["class_%s" % x for x in classes]
    expected = expected_classes + ["id", "total_count"]
    assert out_fields == set(expected)

    for inf, outf in zip(in_features, out_features):
        assert inf["id"] == outf["id"]
        assert inf.geometry().ExportToWkt() == outf.geometry().ExportToWkt()
        total = sum([outf[x] for x in expected_classes])
        assert total == outf["total_count"]

    # Check a couple features
    values = [out_features[3][x] for x in expected_classes]
    assert values == [0, 1, 0, 15, 3, 0]
    values = [out_features[59][x] for x in expected_classes]
    assert values == [151, 3, 0, 1, 18, 0] 
Example 45
Project: bitcoin-arbitrage-trading-bot   Author: mammuth   File: jqvmap.py    MIT License 5 votes vote down vote up
def output_ogr(self, output):
    driver = ogr.GetDriverByName( 'ESRI Shapefile' )
    if os.path.exists( output['file_name'] ):
      driver.DeleteDataSource( output['file_name'] )
    source = driver.CreateDataSource( output['file_name'] )
    layer = source.CreateLayer( self.layer_dfn.GetName(),
                                geom_type = self.layer_dfn.GetGeomType(),
                                srs = self.layer.GetSpatialRef() )

    for field in self.fields:
      fd = ogr.FieldDefn( str(field['name']), field['type'] )
      fd.SetWidth( field['width'] )
      if 'precision' in field:
        fd.SetPrecision( field['precision'] )
      layer.CreateField( fd )

    for geometry in self.geometries:
      if geometry.geom is not None:
        feature = ogr.Feature( feature_def = layer.GetLayerDefn() )
        for index, field in enumerate(self.fields):
          if field['name'] in geometry.properties:
            feature.SetField( index, geometry.properties[field['name']].encode('utf-8') )
          else:
            feature.SetField( index, '' )
        feature.SetGeometryDirectly(
          ogr.CreateGeometryFromWkb(
            shapely.wkb.dumps(
              geometry.geom
            )
          )
        )
        layer.CreateFeature( feature )
        feature.Destroy()

    source.Destroy() 
Example 46
Project: geo-scripts-python   Author: jbranigan   File: detect-utm-zone.py    MIT License 5 votes vote down vote up
def get_bbox(shapefile):
    ''' Gets the bounding box of a shapefile input '''
    driver = ogr.GetDriverByName('ESRI Shapefile')
    data_source = driver.Open(shapefile, 0) # 0 means read, 1 means write
    # Check to see if shapefile is found.
    if data_source is None:
        print 'Could not open %s' % (shapefile)
    else:
        print 'Opened %s' % (shapefile)
        layer = data_source.GetLayer()
        shape_bbox = layer.GetExtent()
        return shape_bbox 
Example 47
Project: qgis-wps4server   Author: 3liz   File: perform_requests.py    GNU General Public License v3.0 5 votes vote down vote up
def testT07ParseExecuteComplexInput(self):
        """Test if Execute with ComplexInput and Output, given directly with input XML request is executed"""
        self._setFromEnv()
        postpywps = pywps.Pywps(pywps.METHOD_POST)
        executeRequestFile = open(os.path.join(
            pywpsPath, "tests", "requests", "wps_execute_request-complexinput-direct.xml"))
        postinputs = postpywps.parseRequest(executeRequestFile)

        postpywps.performRequest(postinputs)
        postxmldom = minidom.parseString(postpywps.response)
        # compare the raster files
        rasterOrig = open(os.path.join(
            pywpsPath, "tests", "datainputs", "dem.tiff"))
        rasterOrigData = rasterOrig.read()
        outputs = postxmldom.getElementsByTagNameNS(self.wpsns, "ComplexData")
        rasterWpsData = base64.decodestring(outputs[1].firstChild.nodeValue)
        self.assertTrue(rasterWpsData, rasterOrigData)

        if os.name != "java":
            # compare the vector files
            gmlDriver = ogr.GetDriverByName("GML")
            origDs = gmlDriver.Open(os.path.join(
                pywpsPath, "tests", "datainputs", "lakes.gml"))

            wpsFile = tempfile.mktemp(prefix="pywps-test")
            wpsFile = open(wpsFile, "w")
            wpsFile.write(postinputs["datainputs"][1]["value"])
            wpsFile.close()
            wpsDs = gmlDriver.Open(wpsFile.name)

            wpslayer = wpsDs.GetLayerByIndex(0)
            origlayer = origDs.GetLayerByIndex(0)

            self.assertTrue(wpslayer.GetFeatureCount(),
                            origlayer.GetFeatureCount())

            # enough  here
            # for f in range(wpslayer.GetFeatureCount()):
            #     origFeature = origlayer.GetFeature(f)
            #     wpsFeature = wpslayer.GetFeature(f)
            #     self.assertTrue(origFeature.Equal(wpsFeature)) 
Example 48
Project: qgis-wps4server   Author: 3liz   File: parser.py    GNU General Public License v3.0 5 votes vote down vote up
def testParseExecuteComplexInputDirectly(self):
        """Test if Execute request is parsed, complex data inputs, given as """

        postpywps = pywps.Pywps(pywps.METHOD_POST)
        executeRequestFile = open(os.path.join(
            pywpsPath, "tests", "requests", "wps_execute_request-complexinput-direct.xml"))
        postinputs = postpywps.parseRequest(executeRequestFile)

        self.assertEquals(postinputs["request"], "execute")
        self.assertTrue("complexprocess" in postinputs["identifier"])
        rasterOrig = open(os.path.join(
            pywpsPath, "tests", "datainputs", "dem.tiff"))
        rasterOrigData = rasterOrig.read()
        rasterWpsData = base64.decodestring(
            postinputs["datainputs"][0]["value"])

        self.assertEquals(rasterOrigData, rasterWpsData)

        gmlDriver = ogr.GetDriverByName("GML")
        origDs = gmlDriver.Open(os.path.join(
            pywpsPath, "tests", "datainputs", "lakes.gml"))

        wpsFile = tempfile.mktemp(prefix="pywps-test")
        wpsFile = open(wpsFile, "w")
        wpsFile.write(postinputs["datainputs"][1]["value"])
        wpsFile.close()
        wpsDs = gmlDriver.Open(wpsFile.name)

        wpslayer = wpsDs.GetLayerByIndex(0)
        origlayer = origDs.GetLayerByIndex(0)

        self.assertTrue(wpslayer.GetFeatureCount(),
                        origlayer.GetFeatureCount())

        # enough  here
        # for f in range(wpslayer.GetFeatureCount()):
        #     origFeature = origlayer.GetFeature(f)
        #     wpsFeature = wpslayer.GetFeature(f)
        #     self.assertTrue(origFeature.Equal(wpsFeature)) 
Example 49
Project: dzetsaka   Author: nkarasiak   File: function_vector.py    GNU General Public License v3.0 5 votes vote down vote up
def saveToShape(self, array, srs, outShapeFile):
        # Parse a delimited text file of volcano data and create a shapefile
        # use a dictionary reader so we can access by field name
        # set up the shapefile driver
        outDriver = ogr.GetDriverByName('ESRI Shapefile')

        # create the data source
        if os.path.exists(outShapeFile):
            outDriver.DeleteDataSource(outShapeFile)
        # Remove output shapefile if it already exists

        # options = ['SPATIALITE=YES'])
        ds = outDriver.CreateDataSource(outShapeFile)

        # create the spatial reference, WGS84

        lyrout = ds.CreateLayer('randomSubset', srs)
        fields = [
            array[1].GetFieldDefnRef(i).GetName() for i in range(
                array[1].GetFieldCount())]

        for f in fields:
            field_name = ogr.FieldDefn(f, ogr.OFTString)
            field_name.SetWidth(24)
            lyrout.CreateField(field_name)

        for k in array:
            lyrout.CreateFeature(k)

        # Save and close the data source
        ds = None 
Example 50
Project: dzetsaka   Author: nkarasiak   File: function_vector.py    GNU General Public License v3.0 5 votes vote down vote up
def saveToShape(array, srs, outShapeFile):
    # Parse a delimited text file of volcano data and create a shapefile
    # use a dictionary reader so we can access by field name
    # set up the shapefile driver
    import ogr
    outDriver = ogr.GetDriverByName('ESRI Shapefile')

    # create the data source
    if os.path.exists(outShapeFile):
        outDriver.DeleteDataSource(outShapeFile)
    # Remove output shapefile if it already exists

    # options = ['SPATIALITE=YES'])
    ds = outDriver.CreateDataSource(outShapeFile)

    # create the spatial reference, WGS84

    lyrout = ds.CreateLayer('randomSubset', srs, ogr.wkbPoint)
    fields = [array[1].GetFieldDefnRef(i).GetName()
              for i in range(array[1].GetFieldCount())]

    for i, f in enumerate(fields):
        field_name = ogr.FieldDefn(f, array[1].GetFieldDefnRef(i).GetType())
        field_name.SetWidth(array[1].GetFieldDefnRef(i).GetWidth())
        lyrout.CreateField(field_name)

    for k in array:
        lyrout.CreateFeature(k)

    # Save and close the data source
    ds = None