Python osgeo.ogr.Open() Examples

The following are code examples for showing how to use osgeo.ogr.Open(). They are extracted from open source Python projects. You can vote up the examples you like or vote down the ones you don't like. You can also save this page to your account.

Example 1
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 7 votes vote down vote up
def convert_pixgwktList_to_wgs84wktList(inputRaster, wktPolygonPixList):
    ## returns # [[GeoWKT, PixWKT], ...]
    wgs84WKTList=[]
    if os.path.isfile(inputRaster):
        srcRaster = gdal.Open(inputRaster)
        targetSR = osr.SpatialReference()
        targetSR.ImportFromWkt(srcRaster.GetProjectionRef())
        geomTransform = srcRaster.GetGeoTransform()

    for wktPolygonPix in wktPolygonPixList:
        geom_wkt_list = pixelWKTToGeoWKT(wktPolygonPix, inputRaster, targetSR='',
                                         geomTransform=geomTransform,
                                         breakMultiPolygonPix=False)

        wgs84WKTList.extend(geom_wkt_list)

    # [[GeoWKT, PixWKT], ...]
    return wgs84WKTList 
Example 2
Project: HistoricalMap   Author: lennepkade   File: function_historical_map.py    (license) View Source Project 6 votes vote down vote up
def polygonize(self,rasterTemp,outShp):
            
            sourceRaster = gdal.Open(rasterTemp)
            band = sourceRaster.GetRasterBand(1)
            driver = ogr.GetDriverByName("ESRI Shapefile")
            # If shapefile already exist, delete it
            if os.path.exists(outShp):
                driver.DeleteDataSource(outShp)
                
            outDatasource = driver.CreateDataSource(outShp)            
            # get proj from raster            
            srs = osr.SpatialReference()
            srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
            # create layer with proj
            outLayer = outDatasource.CreateLayer(outShp,srs)
            # Add class column (1,2...) to shapefile
      
            newField = ogr.FieldDefn('Class', ogr.OFTInteger)
            outLayer.CreateField(newField)
            
            gdal.Polygonize(band, None,outLayer, 0,[],callback=None)  
            
            outDatasource.Destroy()
            sourceRaster=None
            band=None
                    
            
            ioShpFile = ogr.Open(outShp, update = 1)
            
            
            lyr = ioShpFile.GetLayerByIndex(0)
            
            lyr.ResetReading()    

            for i in lyr:
                lyr.SetFeature(i)
            # if area is less than inMinSize or if it isn't forest, remove polygon 
                if i.GetField('Class')!=1:
                    lyr.DeleteFeature(i.GetFID())        
            ioShpFile.Destroy()
                        
            return outShp 
Example 3
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 6 votes vote down vote up
def import_summary_geojson(geojsonfilename, removeNoBuildings=True):
    # driver = ogr.GetDriverByName('geojson')
    datasource = ogr.Open(geojsonfilename, 0)

    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())

    buildingList = []
    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:
            if removeNoBuildings:
                if feature.GetField('BuildingId') != -1:
                    buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
                                  'poly': feature.GetGeometryRef().Clone()})
            else:

                buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
                              'poly': feature.GetGeometryRef().Clone()})

    return buildingList 
Example 4
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 6 votes vote down vote up
def import_chip_geojson(geojsonfilename, ImageId=''):
    # driver = ogr.GetDriverByName('geojson')
    datasource = ogr.Open(geojsonfilename, 0)
    if ImageId=='':
        ImageId = geojsonfilename
    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())

    polys = []
    BuildingId = 0
    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:
            BuildingId = BuildingId + 1
            polys.append({'ImageId': ImageId, 'BuildingId': BuildingId,
                          'poly': feature.GetGeometryRef().Clone()})

    return polys 
Example 5
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 6 votes vote down vote up
def mergePolyList(geojsonfilename):

    multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)
    datasource = ogr.Open(geojsonfilename, 0)

    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())


    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:

            multipolygon.AddGeometry(feature.GetGeometryRef().Clone())

    return multipolygon 
Example 6
Project: rastercube   Author: terrai   File: shputils.py    (license) View Source Project 6 votes vote down vote up
def rasterize_shapefile_like(shpfile, model_raster_fname, dtype, options,
                             nodata_val=0,):
    """
    Given a shapefile, rasterizes it so it has
    the exact same extent as the given model_raster

    `dtype` is a gdal type like gdal.GDT_Byte
    `options` should be a list that will be passed to GDALRasterizeLayers
        papszOptions, like ["ATTRIBUTE=vegetation","ALL_TOUCHED=TRUE"]
    """
    model_dataset = gdal.Open(model_raster_fname)
    shape_dataset = ogr.Open(shpfile)
    shape_layer = shape_dataset.GetLayer()
    mem_drv = gdal.GetDriverByName('MEM')
    mem_raster = mem_drv.Create(
        '',
        model_dataset.RasterXSize,
        model_dataset.RasterYSize,
        1,
        dtype
    )
    mem_raster.SetProjection(model_dataset.GetProjection())
    mem_raster.SetGeoTransform(model_dataset.GetGeoTransform())
    mem_band = mem_raster.GetRasterBand(1)
    mem_band.Fill(nodata_val)
    mem_band.SetNoDataValue(nodata_val)

    # http://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1
    # http://gis.stackexchange.com/questions/31568/gdal-rasterizelayer-doesnt-burn-all-polygons-to-raster
    err = gdal.RasterizeLayer(
        mem_raster,
        [1],
        shape_layer,
        None,
        None,
        [1],
        options
    )
    assert err == gdal.CE_None
    return mem_raster.ReadAsArray() 
Example 7
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 6 votes vote down vote up
def raster2array(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    nodata = band.GetNoDataValue()
    array = band.ReadAsArray()
    
    proj = raster.GetProjection()
    inproj = osr.SpatialReference()
    inproj.ImportFromWkt(proj)
    
    geoTransform = raster.GetGeoTransform()
    minx = geoTransform[0]
    maxy = geoTransform[3]
    maxx = minx + geoTransform[1]*raster.RasterXSize
    miny = maxy + geoTransform[5]*raster.RasterYSize
    extent =  [minx, maxx, miny, maxy]
    pixelSizeXY = [geoTransform[1], geoTransform[5]]
    del raster, band
    return [array, nodata, extent, inproj, pixelSizeXY]

#clip a raster by vector 
Example 8
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 6 votes vote down vote up
def copyproj(src_fn, dst_fn, gt=True):
    """Copy projection and geotransform from one raster file to another
    """
    src_ds = gdal.Open(src_fn, gdal.GA_ReadOnly)
    dst_ds = gdal.Open(dst_fn, gdal.GA_Update)
    dst_ds.SetProjection(src_ds.GetProjection())
    if gt:
        src_gt = np.array(src_ds.GetGeoTransform())
        src_dim = np.array([src_ds.RasterXSize, src_ds.RasterYSize])
        dst_dim = np.array([dst_ds.RasterXSize, dst_ds.RasterYSize])
        #This preserves dst_fn resolution
        if np.any(src_dim != dst_dim):
            res_factor = src_dim/dst_dim.astype(float)
            src_gt[[1, 5]] *= max(res_factor)
            #src_gt[[1, 5]] *= min(res_factor)
            #src_gt[[1, 5]] *= res_factor
        dst_ds.SetGeoTransform(src_gt)
    src_ds = None
    dst_ds = None 
Example 9
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 6 votes vote down vote up
def shp2geom(shp_fn):
    """Extract geometries from input shapefile
    
    Need to handle multi-part geom: http://osgeo-org.1560.x6.nabble.com/Multipart-to-singlepart-td3746767.html
    """
    ds = ogr.Open(shp_fn)
    lyr = ds.GetLayer()
    srs = lyr.GetSpatialRef()
    lyr.ResetReading()
    geom_list = []
    for feat in lyr:
        geom = feat.GetGeometryRef()
        geom.AssignSpatialReference(srs)
        #Duplicate the geometry, or segfault
        #See: http://trac.osgeo.org/gdal/wiki/PythonGotchas
        #g = ogr.CreateGeometryFromWkt(geom.ExportToWkt())
        #g.AssignSpatialReference(srs)
        g = geom_dup(geom)
        geom_list.append(g)
    #geom = ogr.ForceToPolygon(' '.join(geom_list))    
    #Dissolve should convert multipolygon to single polygon 
    #return geom_list[0]
    ds = None
    return geom_list 
Example 10
Project: QGIS-applications   Author: liaduarte   File: GdalTools_utils.py    (license) View Source Project 6 votes vote down vote up
def getVectorFields(vectorFile):
    hds = ogr.Open( unicode(vectorFile).encode('utf8') )
    if hds == None:
      raise UnsupportedOGRFormat()

    fields = []
    names = []

    layer = hds.GetLayer(0)
    defn = layer.GetLayerDefn()

    for i in range(defn.GetFieldCount()):
      fieldDefn = defn.GetFieldDefn(i)
      fieldType = fieldDefn.GetType()
      if fieldType == 0 or fieldType == 2:
        fields.append(fieldDefn)
        names.append(fieldDefn.GetName())

    return (fields, names)

# get raster SRS if possible 
Example 11
Project: dzetsaka   Author: lennepkade   File: mainfunction.py    (license) View Source Project 6 votes vote down vote up
def rasterize(self,inRaster,inShape,inField):
        filename = tempfile.mktemp('.tif')        
        data = gdal.Open(inRaster,gdal.GA_ReadOnly)
        shp = ogr.Open(inShape)
        
        lyr = shp.GetLayer()

        driver = gdal.GetDriverByName('GTiff')
        dst_ds = driver.Create(filename,data.RasterXSize,data.RasterYSize, 1,gdal.GDT_Byte)
        dst_ds.SetGeoTransform(data.GetGeoTransform())
        dst_ds.SetProjection(data.GetProjection())
        OPTIONS = 'ATTRIBUTE='+inField
        gdal.RasterizeLayer(dst_ds, [1], lyr, None,options=[OPTIONS])
        data,dst_ds,shp,lyr=None,None,None,None
        
        
        return filename 
Example 12
Project: beachfront-py   Author: venicegeo   File: vectorize.py    (license) View Source Project 6 votes vote down vote up
def simplify(filename, tolerance=0.00035):
    """ Simplify GeoJSON vector """
    # tolerance is set using GeoJSON map units. Expects decimal degrees, but should work with anything
    if tolerance is None:
        return filename
    logger.info('Updating file %s with simplified geometries' % filename, action='Updating file',
                actee=filename, actor=__name__)
    vs = ogr.Open(filename, 1)  # 1 opens the file in read/write mode, 0 for read-only mode
    layer = vs.GetLayer()
    feat = layer.GetNextFeature()
    while feat is not None:
        geo = feat.geometry()
        simple = geo.Simplify(tolerance)
        feat.SetGeometry(simple)
        layer.SetFeature(feat)
        feat = layer.GetNextFeature()
    layer = None
    vs.Destroy() 
Example 13
Project: beachfront-py   Author: venicegeo   File: test_vectorize.py    (license) View Source Project 6 votes vote down vote up
def test_simplify(self):
        """ Simplify GeoJSON geometries """
        geoimg = download_image(self.test_url)
        geoimg.set_nodata(0)
        lines = vectorize.potrace(geoimg[0] > 9500, close=0)
        fout = os.path.join(os.path.dirname(__file__), 'test.geojson')
        vectorize.save_geojson(lines, fout)

        # check file
        df = ogr.Open(fout)
        layer = df.GetLayer()
        self.assertEqual(layer.GetFeatureCount(), len(lines))
        geom = json.loads(layer.GetNextFeature().ExportToJson())
        self.assertEqual(len(geom['geometry']['coordinates']), len(lines[0]))
        df = None

        # simplify and check file
        vectorize.simplify(fout, tolerance=0.001)
        df = ogr.Open(fout)
        layer = df.GetLayer()
        geom = json.loads(layer.GetNextFeature().ExportToJson())
        self.assertEqual(len(geom['geometry']['coordinates']), 22) 
Example 14
Project: TF-SegNet   Author: mathildor   File: import_data.py    (license) View Source Project 6 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 15
Project: CHaMP_Metrics   Author: SouthForkResearch   File: substrate_raster.py    (license) View Source Project 6 votes vote down vote up
def rasterize_polygons(polygon_shp, template_raster, out_raster, field):
    """ generate a categorical raster based on polygons

    :rtype: None
    :param polygon_shp: input polygon shapefile
    :param template_raster: raster template for cellsize and extent
    :param out_raster: output raster file
    """

    # Source: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

    gdal.UseExceptions()
    # Get template raster specs
    src_ds = gdal.Open(template_raster)
    if src_ds is None:
        print 'Unable to open %s' % template_raster
        sys.exit(1)
    try:
        srcband = src_ds.GetRasterBand(1)
    except RuntimeError, e:
        print 'No band %i found' % 0
        print e
        sys.exit(1)

    # Open the data source and read in the extent
    source_ds = ogr.Open(polygon_shp)
    source_layer = source_ds.GetLayer()

    target_ds = gdal.GetDriverByName('GTiff').Create(out_raster, src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Float32)
    target_ds.SetGeoTransform(src_ds.GetGeoTransform())
    target_ds.SetProjection(src_ds.GetProjection())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(srcband.GetNoDataValue())

    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], source_layer, options=["ATTRIBUTE={}".format(field)]) 
Example 16
Project: HistoricalMap   Author: lennepkade   File: function_historical_map.py    (license) View Source Project 5 votes vote down vote up
def reclassAndFillHole(self,rasterTemp,inClassNumber):
        """ [email protected] Reclass and file hole (with ndimage morphology binary fill hole) for last treatment (need raster and class number)"""
        dst_ds = gdal.Open(rasterTemp,gdal.GA_Update)
        srcband = dst_ds.GetRasterBand(1).ReadAsArray()
        # All data which is not forest is set to 0, so we fill all for the forest only, because it's a binary fill holes.            
        # Set selected class as 1                   
        srcband[srcband != inClassNumber]=0
        srcband[srcband == inClassNumber]=1
        
        srcband = ndimage.morphology.binary_fill_holes(srcband)
        
        dst_ds.GetRasterBand(1).WriteArray(srcband)
        dst_ds.FlushCache()
        return rasterTemp 
Example 17
Project: batch_gps_importer   Author: wondie   File: process_import.py    (license) View Source Project 5 votes vote down vote up
def init_gpx_import(self, gpx_path):
        """
        Initializes the gpx import by setting the gpx path. This method must be
        called outside the class to properly connect signals.
        :param gpx_path: The gpx file path.
        :type gpx_path: String
        """
        # Open GPX file
        if self._file_is_readable(gpx_path):
            self.gpx_file_name = os.path.basename(gpx_path)
            self.file_name = self.gpx_file_name.rstrip('.gpx')
            self.gpx_path = gpx_path
            self.gpx_util = GpxUtil(self.gpx_path)

            data_source = ogr.Open(gpx_path)
            if data_source is not None:
                for feature_type in self.feature_types:
                    if STOP_IMPORT:
                        return
                    self.feature_type = feature_type.encode('utf-8')
                    self.ogr_layer = data_source.GetLayerByName(
                        self.feature_type
                    )

                    if self.ogr_layer.GetFeatureCount() > 0:
                        try:
                            self.gpx_to_point_list()
                        except Exception as ex:
                            self.error_type = ex
                            self.add_progress()
                        if STOP_IMPORT:
                            return
                        self.save_valid_folders()
                        self.save_invalid_folders() 
Example 18
Project: Planet-GEE-Pipeline-GUI   Author: samapriya   File: ogr2ft.py    (license) View Source Project 5 votes vote down vote up
def _get_ft_ds():
    refresh_token = OAuth2().get_refresh_token()
    ft_driver = ogr.GetDriverByName('GFT')

    return ft_driver.Open('GFT:refresh=' + refresh_token, True) 
Example 19
Project: utilities   Author: SpaceNetChallenge   File: externalVectorProcessing.py    (license) View Source Project 5 votes vote down vote up
def buildTindex(rasterFolder, rasterExtention='.tif'):
    rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
    print(rasterList)

    print(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))

    memDriver = ogr.GetDriverByName('MEMORY')
    gTindex = memDriver.CreateDataSource('gTindex')
    srcImage = gdal.Open(rasterList[0])
    spat_ref = osr.SpatialReference()
    spat_ref.SetProjection(srcImage.GetProjection())
    gTindexLayer = gTindex.CreateLayer("gtindexlayer", spat_ref, geom_type=ogr.wkbPolygon)

    # Add an ID field
    idField = ogr.FieldDefn("location", ogr.OFTString)
    gTindexLayer.CreateField(idField)

    # Create the feature and set values
    featureDefn = gTindexLayer.GetLayerDefn()



    for rasterFile in rasterList:
        srcImage = gdal.Open(rasterFile)

        geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)

        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(polyToCut)
        feature.SetField("location", rasterFile)
        gTindexLayer.CreateFeature(feature)
        feature = None


    return gTindex, gTindexLayer 
Example 20
Project: utilities   Author: SpaceNetChallenge   File: externalVectorProcessing.py    (license) View Source Project 5 votes vote down vote up
def createTiledGeoJsonFromSrc(rasterFolderLocation, vectorSrcFile, geoJsonOutputDirectory, rasterTileIndex='',
                              geoJsonPrefix='GEO', rasterFileExtenstion='.tif',
                              rasterPrefixToReplace='PAN'
                              ):
    if rasterTileIndex == '':
        gTindex, gTindexLayer = buildTindex(rasterFolderLocation, rasterExtention=rasterFileExtenstion)
    else:
        gTindex = ogr.Open(rasterTileIndex,0)
        gTindexLayer = gTindex.GetLayer()

    shapeSrc = ogr.Open(vectorSrcFile,0)
    chipSummaryList = []
    for feature in gTindexLayer:
        featureGeom = feature.GetGeometryRef()
        rasterFileName = feature.GetField('location')
        rasterFileBaseName = os.path.basename(rasterFileName)
        outGeoJson = rasterFileBaseName.replace(rasterPrefixToReplace, geoJsonPrefix)
        outGeoJson = outGeoJson.replace(rasterFileExtenstion, '.geojson')
        outGeoJson = os.path.join(geoJsonOutputDirectory, outGeoJson)

        gT.clipShapeFile(shapeSrc, outGeoJson, featureGeom, minpartialPerc=0.0, debug=False)
        imageId = rasterFileBaseName.replace(rasterPrefixToReplace+"_", "")
        chipSummary = {'chipName': rasterFileName,
                           'geoVectorName': outGeoJson,
                           'imageId': os.path.splitext(imageId)[0]}

        chipSummaryList.append(chipSummary)

    return chipSummaryList 
Example 21
Project: utilities   Author: SpaceNetChallenge   File: labelTools.py    (license) View Source Project 5 votes vote down vote up
def createTruthPixelLinePickle(truthLineFile, pickleLocation=''):
    if pickleLocation=='':
        extension = os.path.splitext(truthLineFile)[1]
        pickleLocation = truthLineFile.replace(extension, 'Pixline.p')
    if truthLineFile != '':
        # get Source Line File Information
        shapef = ogr.Open(truthLineFile, 0)
        truthLayer = shapef.GetLayer()
        pt1X = []
        pt1Y = []
        pt2X = []
        pt2Y = []
        for tmpFeature in truthLayer:
            tmpGeom = tmpFeature.GetGeometryRef()
            for i in range(0, tmpGeom.GetPointCount()):
                pt = tmpGeom.GetPoint(i)

                if i == 0:
                    pt1X.append(pt[0])
                    pt1Y.append(pt[1])
                elif i == 1:
                    pt2X.append(pt[0])
                    pt2Y.append(pt[1])

        lineData = {'pt1X': np.asarray(pt1X),
                    'pt1Y': np.asarray(pt1Y),
                    'pt2X': np.asarray(pt2X),
                    'pt2Y': np.asarray(pt2Y)
                    }

        with open(pickleLocation, 'wb') as f:
            pickle.dump(lineData, f)
            # get Source Line File Information 
Example 22
Project: utilities   Author: SpaceNetChallenge   File: labelTools.py    (license) View Source Project 5 votes vote down vote up
def createTruthPixelPolyPickle(truthPoly, pickleLocation=''):
    # returns dictionary with list of minX, maxX, minY, maxY

    if pickleLocation=='':
        extension = os.path.splitext(truthPoly)[1]
        pickleLocation = truthPoly.replace(extension, 'PixPoly.p')
    if truthPoly != '':
        # get Source Line File Information
        shapef = ogr.Open(truthPoly, 0)
        truthLayer = shapef.GetLayer()
        envList = []

        for tmpFeature in truthLayer:
            tmpGeom = tmpFeature.GetGeometryRef()
            env = tmpGeom.GetEvnelope()
            envList.append(env)

        envArray = np.asarray(envList)
        envelopeData = {'minX': envArray[:,0],
                        'maxX': envArray[:,1],
                        'minY': envArray[:,2],
                        'maxY': envArray[:,3]
                        }


        with open(pickleLocation, 'wb') as f:
            pickle.dump(envelopeData, f)
            # get Source Line File Information 
Example 23
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 5 votes vote down vote up
def clipRaster(raster, newRaster, vector):
    raster = gdal.Open(raster)
    
    vect = ogr.Open(vector)
    lyr = vect.GetLayer()
    ext = lyr.GetExtent()
    
    gTrans = raster.GetGeoTransform()
    #get the x start of the left most pixel
    xlP = int((ext[0] - gTrans[0])/gTrans[1])*gTrans[1] - abs(gTrans[0])
    #get the x end of the right most pixel
    xrP = math.ceil((ext[1] - gTrans[0])/gTrans[1])*gTrans[1] - abs(gTrans[0])
    #get the y start of the upper most pixel
    yuP = abs(gTrans[3]) - int((gTrans[3] - ext[3])/gTrans[5])*gTrans[5]
    #get the y end of the lower most pixel
    ylP = abs(gTrans[3]) - math.floor((gTrans[3] - ext[2])/gTrans[5])*gTrans[5]
        
    gdal.Translate('tmp.tif', raster, projWin = [xlP, yuP, xrP, ylP])
    del raster
    tRas = gdal.Open('tmp.tif')
    band = tRas.GetRasterBand(1)
    noDat = band.GetNoDataValue()
    # store a copy before rasterize
    fullRas = band.ReadAsArray().astype(numpy.float)
    
    gdal.RasterizeLayer(tRas, [1], lyr, None, None, [-9999], ['ALL_TOUCHED=TRUE']) # now tRas is modified
    
    finRas = tRas.GetRasterBand(1).ReadAsArray().astype(numpy.float)

    for i, row in enumerate(finRas):
        for j, col in enumerate(row):
            if col == -9999.:
                finRas[i, j] = fullRas[i, j]
            else:
                finRas[i, j] = noDat

    array2raster(newRaster, 'tmp.tif', finRas, noDat, "float32")
    os.remove('tmp.tif')
    del fullRas, finRas, band, tRas

# create a reference raster with random values 
Example 24
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 5 votes vote down vote up
def array2raster(newRaster, RefRaster, array, noData, dataType):
    #data type conversion
    NP2GDAL_CONVERSION = { "uint8": 1, "int8": 1, "uint16": 2, "int16": 3, 
                          "uint32": 4, "int32": 5, "float32": 6, "float64": 7,
                          "complex64": 10, "complex128": 11,
                         }
    #get info from reference raster
    rfRaster = gdal.Open(RefRaster)
    geotransform = rfRaster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = array.shape[1]
    rows = array.shape[0]
    #create new raster
    outRaster = gdal.GetDriverByName('GTiff').Create(newRaster, cols, rows,1, NP2GDAL_CONVERSION[dataType])
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    #write array to band
    outband = outRaster.GetRasterBand(1)
    outband.SetNoDataValue(noData)
    outband.WriteArray(array)
    #define new raster projection
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(rfRaster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    #write raster
    outband.FlushCache()
    del rfRaster 
Example 25
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def shp_dict(shp_fn, fields=None, geom=True):
    """Get a dictionary for all features in a shapefile
    
    Optionally, specify fields
    """
    from pygeotools.lib import timelib
    ds = ogr.Open(shp_fn)
    lyr = ds.GetLayer()
    nfeat = lyr.GetFeatureCount()
    print('%i input features\n' % nfeat)
    if fields is None:
        fields = shp_fieldnames(lyr)
    d_list = []
    for n,feat in enumerate(lyr):
        d = {}
        if geom:
            geom = feat.GetGeometryRef()
            d['geom'] = geom
        for f_name in fields:
            i = str(feat.GetField(f_name))
            if 'date' in f_name:
                # date_f = f_name
                #If d is float, clear off decimal
                i = i.rsplit('.')[0]
                i = timelib.strptime_fuzzy(str(i))
            d[f_name] = i
        d_list.append(d)
    #d_list_sort = sorted(d_list, key=lambda k: k[date_f])
    return d_list 
Example 26
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def clip_raster_by_shp(dem_fn, shp_fn):
    import subprocess
    #This is ok when writing to outdir, but clip_raster_by_shp.sh writes to raster dir
    #try:
    #    with open(dem_fn) as f: pass
    #except IOError as e:
    cmd = ['clip_raster_by_shp.sh', dem_fn, shp_fn]
    print(cmd)
    subprocess.call(cmd, shell=False)
    dem_clip_fn = os.path.splitext(dem_fn)[0]+'_shpclip.tif'
    dem_clip_ds = gdal.Open(dem_clip_fn, gdal.GA_ReadOnly)
    return dem_clip_ds

#Hack
#extent is xmin ymin xmax ymax 
Example 27
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def wgs84_to_egm96(dem_ds, geoid_dir=None):
    from pygeotools.lib import warplib

    #Check input dem_ds - make sure WGS84

    geoid_dir = os.getenv('ASP_DATA')
    if geoid_dir is None:
        print("No geoid directory available")
        print("Set ASP_DATA or specify")
    
    egm96_fn = geoid_dir+'/geoids-1.1/egm96-5.tif' 
    try:
        open(egm96_fn)
    except IOError:
        sys.exit("Unable to find "+egm96_fn)
    egm96_ds = gdal.Open(egm96_fn)

    #Warp egm96 to match the input ma
    ds_list = warplib.memwarp_multi([dem_ds, egm96_ds], res='first', extent='first', t_srs='first') 

    #Try two-step with extent/res in wgs84
    #ds_list = warplib.memwarp_multi([dem_ds, egm96_ds], res='first', extent='intersection', t_srs='last') 
    #ds_list = warplib.memwarp_multi([dem_ds, ds_list[1]], res='first', extent='first', t_srs='first')

    print("Extracting ma from dem and egm96 ds")
    from pygeotools.lib import iolib
    dem = iolib.ds_getma(ds_list[0])
    egm96 = iolib.ds_getma(ds_list[1])

    print("Removing offset")
    dem_egm96 = dem - egm96
   
    return dem_egm96

#Run ASP dem_geoid adjustment utility
#Note: this is multithreaded 
Example 28
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def dem_geoid(dem_fn):
    out_prefix = os.path.splitext(dem_fn)[0]
    adj_fn = out_prefix +'-adj.tif'
    if not os.path.exists(adj_fn):
        import subprocess
        cmd_args = ["-o", out_prefix, dem_fn]
        cmd = ['dem_geoid'] + cmd_args
        #cmd = 'dem_geoid -o %s %s' % (out_prefix, dem_fn)
        print(' '.join(cmd))
        subprocess.call(cmd, shell=False)
    adj_ds = gdal.Open(adj_fn, gdal.GA_ReadOnly)
    #from pygeotools.lib import iolib
    #return iolib.ds_getma(adj_ds, 1)
    return adj_ds 
Example 29
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 5 votes vote down vote up
def dem_geoid_offsetgrid(dem_fn):
    ds = gdal.Open(dem_fn)
    out_fn = os.path.splitext(dem_fn)[0]+'_EGM2008offset.tif'
    o = dem_geoid_offsetgrid_ds(ds, out_fn)
    return o

#Note: funcitonality with masking needs work 
Example 30
Project: pypgroutingloader   Author: danieluct   File: pgrouting_distance_isobands.py    (license) View Source Project 5 votes vote down vote up
def __init__(self, params, zField='z'):
        ds = ogr.Open(connString)
        layer = ds.ExecuteSQL(OGR_SQL.format(*params))
        extent = layer.GetExtent()

        
        self.proj = layer.GetSpatialRef()
        self.geotransform = []        
        self.x = []
        self.y = []
        self.vals = []

        xMin, xMax, yMin, yMax = extent
        xSize, ySize = abs(xMax - xMin) / 0.0003, abs(yMin - yMax) / 0.0003
        self.size = xSize, ySize
        
        self.geotransform = [xMin, (xMax - xMin) / xSize, 0,
                             yMax, 0, (yMin - yMax) / ySize]

        
        feature = layer.GetNextFeature()        
        if feature.GetFieldIndex(zField) == -1:
            raise Exception('zField is not valid: ' + zField)
        while feature:  
            geometry = feature.GetGeometryRef()
            self.x.append(geometry.GetX())
            self.y.append(geometry.GetY()) 
            self.vals.append(feature.GetField(zField))  
            feature = layer.GetNextFeature()
        ds.Destroy() 
Example 31
Project: layman   Author: CCSS-CZ   File: __init__.py    (license) View Source Project 5 votes vote down vote up
def get_gis_attributes(self,fileName, attrs):
        """Append more gis attributes of given file to attrs hash
        """
        logging.debug("[FileMan][get_gis_attributes] Params: fileName: %s, attrs: %s" % (fileName, repr(attrs)) )       
 
        # try vector
        ds = ogr.Open(fileName)
       
        # opening vector success
        if ds:
            logging.debug("[FileMan][get_gis_attributes] ogr.Open() O.K" )
            attrs = self._get_vector_attributes(ds,attrs)
            logging.debug("[FileMan][get_gis_attributes] Identified VECTOR attributes: %s" % repr(attrs) )       

        # try raster
        else:
            logging.debug("[FileMan][get_gis_attributes] ogr.Open() Failed" )
            ds = gdal.Open(fileName)

            # opening raster success
            if ds:
                logging.debug("[FileMan][get_gis_attributes] gdal.Open() O.K." )
                attrs = self._get_raster_attributes( ds,attrs)
                logging.debug("[FileMan][get_gis_attributes] Identified RASTER attributes: %s" % repr(attrs) )       
            # no gis data
            else:
                logging.debug("[FileMan][get_gis_attributes] gdal.Open() Failed" )
                logging.debug("[FileMan][get_gis_attributes] No attributes identified for file %s" % fileName )       

        return attrs 
Example 32
Project: layman   Author: CCSS-CZ   File: __init__.py    (license) View Source Project 5 votes vote down vote up
def updateData(self, dataName, workspace, fsUserDir, fsGroupDir,
                   dbSchema, fileName):
        """Update data - database or file system -
           from new shape or raster file
        """

        filePath = os.path.realpath(os.path.join(fsUserDir, fileName))

        from osgeo import ogr
        ds = ogr.Open(filePath)
        data_type = None

        # VECTOR
        if ds:

            # Import to DB
            from layman.layed.dbman import DbMan
            dbm = DbMan(self.config)
            dbm.updateVectorFile(filePath, dbSchema, dataName)
            data_type = "vector"

        # RASTER
        else:
            from osgeo import gdal
            ds = gdal.Open(filePath)
            if ds:
                self.updateRasterFile(workspace, filePath)
                data_type = "raster"
                return

        if not data_type:
            raise LaymanError(500,
                              "Data type (raster or vector) not recognized")

    ### STYLES ### 
Example 33
Project: wa   Author: wateraccounting   File: raster_conversions.py    (license) View Source Project 5 votes vote down vote up
def Open_array_info(filename=''):

    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 34
Project: wa   Author: wateraccounting   File: raster_conversions.py    (license) View Source Project 5 votes vote down vote up
def Open_tiff_array(filename='', band=''):	

    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 35
Project: wa   Author: wateraccounting   File: raster_conversions.py    (license) View Source Project 5 votes vote down vote up
def clip_data(input_file, latlim, lonlim):
    """
    Clip the data to the defined extend of the user (latlim, lonlim) or to the
    extend of the DEM tile

    Keyword Arguments:
    input_file -- output data, output of the clipped dataset
    latlim -- [ymin, ymax]
    lonlim -- [xmin, xmax]
    """
    try:			
        if input_file.split('.')[-1] == 'tif':
            dest_in = gdal.Open(input_file)               					
        else:
            dest_in = input_file
    except:
        dest_in = input_file
    
    # Open Array
    data_in = dest_in.GetRasterBand(1).ReadAsArray()	

    # Define the array that must remain
    Geo_in = dest_in.GetGeoTransform()
    Geo_in = list(Geo_in)			
    Start_x = np.max([int(np.ceil(((lonlim[0]) - Geo_in[0])/ Geo_in[1])),0])   				
    End_x = np.min([int(np.floor(((lonlim[1]) - Geo_in[0])/ Geo_in[1])),int(dest_in.RasterXSize)])				
				
    Start_y = np.max([int(np.floor((Geo_in[3] - latlim[1])/ -Geo_in[5])),0])
    End_y = np.min([int(np.ceil(((latlim[0]) - Geo_in[3])/Geo_in[5])), int(dest_in.RasterYSize)])	

    #Create new GeoTransform
    Geo_in[0] = Geo_in[0] + Start_x * Geo_in[1]
    Geo_in[3] = Geo_in[3] + Start_y * Geo_in[5]
    Geo_out = tuple(Geo_in)
				
    data = np.zeros([End_y - Start_y, End_x - Start_x])				

    data = data_in[Start_y:End_y,Start_x:End_x] 
    dest_in = None
    
    return(data, Geo_out) 
Example 36
Project: dzetsaka   Author: lennepkade   File: function_vector.py    (license) View Source Project 5 votes vote down vote up
def __init__(self,inShape,inField,outValidation,outTrain,number=50,percent=True):
        """
        inShape : str path file (e.g. '/doc/ref.shp')
        inField : string column name (e.g. 'class')
        outValidation : str path of shp output file (e.g. '/tmp/valid.shp')
        outTrain : str path of shp output file (e.g. '/tmp/train.shp')
        """
        if percent:
            number = number / 100.0
        else:
            number = int(number)
            
        lyr = ogr.Open(inShape)
        lyr1 = lyr.GetLayer()
        FIDs= sp.zeros(lyr1.GetFeatureCount(),dtype=int)
        Features = []
        #unselFeat = []
        #current = 0
        
        for i,j in enumerate(lyr1):
            #print j.GetField(inField)
            FIDs[i] = j.GetField(inField)
            Features.append(j)
            #current += 1
        srs = lyr1.GetSpatialRef()
        lyr1.ResetReading()
        
        ## 
        if percent:
            validation,train = train_test_split(Features,test_size=number,train_size=1-number,stratify=FIDs)
        else:
            validation,train = train_test_split(Features,test_size=number,stratify=FIDs)
        
        self.saveToShape(validation,srs,outValidation)
        self.saveToShape(train,srs,outTrain) 
Example 37
Project: dzetsaka   Author: lennepkade   File: dzetsaka.py    (license) View Source Project 5 votes vote down vote up
def saveConfusion(self):
        """
        Open window and save confusion shown in qtableview
        """
        fileName = QFileDialog.getSaveFileName(self.dockwidget, "Select output file",self.lastSaveDir,"CSV (*.csv)")
        self.rememberLastSaveDir(fileName)
        fileName,fileExtension=os.path.splitext(fileName)
        if fileExtension != '.csv':
            fileName=fileName+'.csv'
     
        # save to CSV
        try :
            sp.savetxt(fileName,self.lastConfusionMatrix ,delimiter=',',fmt='%1.4d')
        except : 
            QtGui.QMessageBox.warning(self, 'Missing confusion matrix ? ', 'Cannot save confusion matrix. Are you sure to have generated it before ?', QtGui.QMessageBox.Ok) 
Example 38
Project: beachfront-py   Author: venicegeo   File: mask.py    (license) View Source Project 5 votes vote down vote up
def open_vector(filename, layer=''):
    """ Fetch wfs features within bounding box """
    logger.info('Opening %s as vector file' % filename, action='Open file', actee=filename, actor=__name__)
    ds = ogr.Open(filename)
    if layer == '':
        layer = ds.GetLayer(0)
    else:
        layer = ds.GetLayer(layer)
    return (ds, layer) 
Example 39
Project: TF-SegNet   Author: mathildor   File: import_data.py    (license) View Source Project 5 votes vote down vote up
def importBygg(to_cur, filename):
    print("Importing:", filename)
    dataSource = ogr.Open(filename)
    dataLayer = dataSource.GetLayer(0)
    print("dataLayer")
    print(dataLayer)

    #Genererer id'er fortloopende
    for id, feature in enumerate(dataLayer):
        geom = feature.GetGeometryRef()
        if not geom:
            continue
        geom.FlattenTo2D()
        print("Feature")
        print(feature)

        for i in range(1, feature.GetFieldCount()):
            field = feature.GetDefnRef().GetFieldDefn(i).GetName()
            if( i == 4):
                continue
            #print(field.encode('utf-8'))

        byggtyp = feature.GetField("BYGGTYP_NB")
        #poly_id = feature.GetField("LOKALID   ")
        objtype = feature.GetField("OBJTYPE")

        to_tuple = (id, objtype, byggtyp, 'SRID=25832;' + geom.ExportToWkt())

        to_cur.execute("""INSERT into ar_bygg (id, objtype,  byggtyp, geom)
                        SELECT %s, %s, %s, ST_GeometryFromText(%s);""",
                    to_tuple)

    to_conn.commit()
    dataSource.Destroy()

# to_conn = psycopg2.connect("host=localhost port=5433 dbname=ar-bygg-ostfold user=postgres password=24Pils") 
Example 40
Project: HistoricalMap   Author: lennepkade   File: function_historical_map.py    (license) View Source Project 4 votes vote down vote up
def postClassRaster(self,inRaster,sieveSize,inClassNumber,outShp):        
        """ [email protected] Sieve size with gdal.Sieve() fiunction, them reclass to delete unwanted labels """
        
        try:
            rasterTemp = tempfile.mktemp('.tif')
                    
            datasrc = gdal.Open(inRaster)
            srcband = datasrc.GetRasterBand(1)
            data,im=dataraster.open_data_band(inRaster)        
            
            drv = gdal.GetDriverByName('GTiff')
            dst_ds = drv.Create(rasterTemp,datasrc.RasterXSize,datasrc.RasterXSize,1,gdal.GDT_Byte)
            
            dst_ds.SetGeoTransform(datasrc.GetGeoTransform())
            dst_ds.SetProjection(datasrc.GetProjection())
        
            dstband=dst_ds.GetRasterBand(1)
            
            
            
            def sieve(srcband,dstband,sieveSize):
                gdal.SieveFilter(srcband,None,dstband,sieveSize,8)
            
            pixelSize = datasrc.GetGeoTransform()[1] #get pixel size
            pixelSieve = int(sieveSize/(pixelSize*pixelSize)) #get number of pixel to sieve
            
            sieve(srcband,dstband,pixelSieve)
            
            dst_ds = None # close destination band
            
            self.Progress.addStep()
            

            rasterTemp = self.reclassAndFillHole(rasterTemp,inClassNumber)
            
            
            self.Progress.addStep()
            
            
            
        except:
            QgsMessageLog.logMessage("Cannot sieve with raster function")

        
        try :
            outShp = self.polygonize(rasterTemp,outShp) # vectorize raster
        except :
            self.Progress.reset()    
        
        self.Progress.addStep()        
        
        self.Progress.reset()
        return outShp 
Example 41
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 4 votes vote down vote up
def pixelToGeoCoord(xPix, yPix, inputRaster, sourceSR='', geomTransform='', targetSR=''):
    # If you want to garuntee lon lat output, specify TargetSR  otherwise, geocoords will be in image geo reference
    # targetSR = osr.SpatialReference()
    # targetSR.ImportFromEPSG(4326)
    # Transform can be performed at the polygon level instead of pixel level

    if targetSR =='':
        performReprojection=False
        targetSR = osr.SpatialReference()
        targetSR.ImportFromEPSG(4326)
    else:
        performReprojection=True

    if geomTransform=='':
        srcRaster = gdal.Open(inputRaster)
        geomTransform = srcRaster.GetGeoTransform()

        source_sr = osr.SpatialReference()
        source_sr.ImportFromWkt(srcRaster.GetProjectionRef())




    geom = ogr.Geometry(ogr.wkbPoint)
    xOrigin = geomTransform[0]
    yOrigin = geomTransform[3]
    pixelWidth = geomTransform[1]
    pixelHeight = geomTransform[5]

    xCoord = (xPix * pixelWidth) + xOrigin
    yCoord = (yPix * pixelHeight) + yOrigin
    geom.AddPoint(xCoord, yCoord)


    if performReprojection:
        if sourceSR=='':
            srcRaster = gdal.Open(inputRaster)
            sourceSR = osr.SpatialReference()
            sourceSR.ImportFromWkt(srcRaster.GetProjectionRef())
        coord_trans = osr.CoordinateTransformation(sourceSR, targetSR)
        geom.Transform(coord_trans)

    return (geom.GetX(), geom.GetY()) 
Example 42
Project: utilities   Author: SpaceNetChallenge   File: labelTools.py    (license) View Source Project 4 votes vote down vote up
def convertLabelStringToPoly(shapeFileSrc, outGeoJSon, labelType='Airplane'):

        shapeSrc = ogr.Open(shapeFileSrc)
        source_layer = shapeSrc.GetLayer()
        source_srs = source_layer.GetSpatialRef()
        # Create the output Layer
        outDriver = ogr.GetDriverByName("geojson")
        if os.path.exists(outGeoJSon):
            outDriver.DeleteDataSource(outGeoJSon)


        outDataSource = outDriver.CreateDataSource(outGeoJSon)
        outLayer = outDataSource.CreateLayer("groundTruth", source_srs, geom_type=ogr.wkbPolygon)
        # Add input Layer Fields to the output Layer
        inLayerDefn = source_layer.GetLayerDefn()
        for i in range(0, inLayerDefn.GetFieldCount()):
            fieldDefn = inLayerDefn.GetFieldDefn(i)
            outLayer.CreateField(fieldDefn)
        outLayer.CreateField(ogr.FieldDefn("Length_m", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("Width_m", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("Aspect(L/W)", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("compassDeg", ogr.OFTReal))

        outLayerDefn = outLayer.GetLayerDefn()
        for inFeature in source_layer:

            outFeature = ogr.Feature(outLayerDefn)

            for i in range(0, inLayerDefn.GetFieldCount()):
                outFeature.SetField(inLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))

            geom = inFeature.GetGeometryRef()
            if labelType == 'Airplane':
                poly, Length, Width, Aspect, Direction = evaluateLineStringPlane(geom, label='Airplane')
            elif labelType == 'Boat':
                poly, Length, Width, Aspect, Direction = evaluateLineStringBoat(geom, label='Boat')

            outFeature.SetGeometry(poly)
            outFeature.SetField("Length_m", Length)
            outFeature.SetField("Width_m", Width)
            outFeature.SetField("Aspect(L/W)", Aspect)
            outFeature.SetField("compassDeg", Direction)

            outLayer.CreateFeature(outFeature) 
Example 43
Project: utilities   Author: SpaceNetChallenge   File: labelTools.py    (license) View Source Project 4 votes vote down vote up
def geoJsonToSBD(annotationName_cls, annotationName_inst, geoJson, rasterSource):

    #Print raster file name
    my_raster_source = rasterSource
    print("Raster directory : ",my_raster_source)
    srcRaster = gdal.Open(my_raster_source)

    
    my_vector_source = geoJson
    print("Vector directory : ",my_vector_source)

    
    #Call main functions to create label datafor cls
    my_cls_segmentation = createClassSegmentation(my_raster_source, my_vector_source, npDistFileName='', units='pixels')
    my_cls_boundaries =  createClassBoundaries(my_raster_source, my_vector_source, npDistFileName='', units='pixels')
    my_cls_categories = createClassCategoriesPresent(my_vector_source)

    #Call main functions to create label datafor inst
    my_inst_segmentation = createInstanceSegmentation(my_raster_source, my_vector_source)
    my_inst_boundaries = createInstanceBoundaries(my_raster_source, my_vector_source)
    my_inst_categories = createInstanceCategories(my_vector_source)

    #Wraps for cls struct
    cls_boundaries_wrap = np.array([my_cls_boundaries])
    cls_categories_wrap = my_cls_categories

    #Wraps for inst struct
    inst_boundaries_wrap = np.array([my_inst_boundaries])
    inst_categories_wrap = my_inst_categories

    #Create a class struct
    GTcls = {'Segmentation': my_cls_segmentation , 'Boundaries': cls_boundaries_wrap, 'CategoriesPresent': cls_categories_wrap}


    #Create the instance struct
    GTinst = {'Segmentation': my_inst_segmentation , 'Boundaries': inst_boundaries_wrap, 'Categories': inst_categories_wrap}

    #Save the files
    scipy.io.savemat(annotationName_cls,{'GTcls': GTcls})
    scipy.io.savemat(annotationName_inst,{'GTinst': GTinst})

    print("Done with "+str())

    entry = {'rasterFileName': my_raster_source,
             'geoJsonFileName': geoJson,
             'annotationName' : annotationName_cls,
             'annotationName_cls': annotationName_cls,
             'annotationName_inst':annotationName_inst,
             'width': srcRaster.RasterXSize,
             'height': srcRaster.RasterYSize,
             'depth' : srcRaster.RasterCount,
             'basename': os.path.splitext(os.path.basename(my_raster_source))[0]
             }

    return entry 
Example 44
Project: rastercube   Author: terrai   File: regions.py    (license) View Source Project 4 votes vote down vote up
def __init__(self, fname):
        ds = ogr.Open(fname)
        assert ds is not None, "Couldn't open %s" % fname
        layer = ds.GetLayer()
        assert layer is not None

        WGS84_SR = osr.SpatialReference()
        WGS84_SR.ImportFromEPSG(4326)

        shape_sr = osr.SpatialReference()
        shape_sr.ImportFromWkt(layer.GetSpatialRef().ExportToWkt())

        transform = osr.CoordinateTransformation(shape_sr, WGS84_SR)

        # http://geoinformaticstutorial.blogspot.ch/2012/10/accessing-vertices-from-polygon-with.html
        polygons = {}
        for feature in layer:
            attr = feature.items()
            newattr = {}
            # some attributes contains unicode character. ASCIIfy everything
            # TODO: A bit brutal, but easy...
            for k, v in attr.items():
                newk = k.decode('ascii', errors='ignore')
                newattr[newk] = v
            attr = newattr

            geometry = feature.GetGeometryRef()
            assert geometry.GetGeometryName() == 'POLYGON'
            # A ring is a polygon in shapefiles
            ring = geometry.GetGeometryRef(0)
            assert ring.GetGeometryName() == 'LINEARRING'
            # The ring duplicates the last point, so for the polygon to be
            # closed, last point should equal first point
            npoints = ring.GetPointCount()
            points = [ring.GetPoint(i) for i in xrange(npoints)]
            points = [transform.TransformPoint(*p) for p in points]
            # third column is elevation - discard
            points = np.array(points)[:, :2]
            # swap (lng, lat) to (lat, lng)
            points = points[:, ::-1]
            assert np.allclose(points[-1], points[0])

            name = attr['name']
            polygons[name] = points

        self.polygons = polygons 
Example 45
Project: rastercube   Author: terrai   File: shputils.py    (license) View Source Project 4 votes vote down vote up
def load_polygons_from_shapefile(filename, target_sr):
    """
    Loads the given shapefiles, reprojects the polygons it contains in the
    given target spatialreference.
    Returns:
        polygons: A list of polygons (list of points)
        attributes: A list of attributes (dict of strings)
    """
    shape = ogr.Open(filename)
    assert shape, "Couldn't open %s" % filename
    assert shape.GetLayerCount() == 1
    layer = shape.GetLayer()
    nfeatures = layer.GetFeatureCount()

    shape_sr = osr.SpatialReference()
    shape_sr.ImportFromWkt(layer.GetSpatialRef().ExportToWkt())

    # Transform from shape to image coordinates
    transform = osr.CoordinateTransformation(shape_sr, target_sr)

    # http://geoinformaticstutorial.blogspot.ch/2012/10/accessing-vertices-from-polygon-with.html
    polygons = []
    attributes = []
    for i in xrange(nfeatures):
        feature = layer.GetFeature(i)
        attr = feature.items()
        newattr = {}
        # some attributes contains unicode character. ASCIIfy everything
        # TODO: A bit brutal, but easy...
        for k, v in attr.items():
            newk = k.decode('ascii', errors='ignore')
            newattr[newk] = v
        attr = newattr

        geometry = feature.GetGeometryRef()
        assert geometry.GetGeometryName() == 'POLYGON'
        # A ring is a polygon in shapefiles
        ring = geometry.GetGeometryRef(0)
        assert ring.GetGeometryName() == 'LINEARRING'
        # The ring duplicates the last point, so for the polygon to be closed,
        # last point should equal first point
        npoints = ring.GetPointCount()
        points = [ring.GetPoint(i) for i in xrange(npoints)]
        points = [transform.TransformPoint(*p) for p in points]
        points = np.array(points)[:, :2]  # third column is elevation - discard
        # swap (lng, lat) to (lat, lng)
        points = points[:, ::-1]
        assert np.allclose(points[-1], points[0])
        polygons.append(points)
        attributes.append(attr)

    return polygons, attributes 
Example 46
Project: chorospy   Author: spyrostheodoridis   File: bioFunc.py    (license) View Source Project 4 votes vote down vote up
def makeDensityRaster(speciesOcc, inVector, pixelSize, outRas, noData):
    srcVector = ogr.Open(inVector)
    srcLayer = srcVector.GetLayer()
    srs = srcLayer.GetSpatialRef()
    # if the layer is not wgs84
    if srs.GetAttrValue("AUTHORITY", 1) != '4326':
        print('Layer projection should be WGS84!')
        return

    xMin, xMax, yMin, yMax = srcLayer.GetExtent()

    # Create the destination data source
    xRes = int((xMax - xMin) / pixelSize)
    yRes = int((yMax - yMin) / pixelSize)
    targetRas = gdal.GetDriverByName('GTiff').Create(outRas, xRes, yRes, 1, 6) # 6 == float
    targetRas.SetGeoTransform((xMin, pixelSize, 0, yMax, 0, -pixelSize))
    band = targetRas.GetRasterBand(1)
    band.SetNoDataValue(noData)

    # Rasterize clips the raster band
    gdal.RasterizeLayer(targetRas, [1], srcLayer, None, None, [0], ['ALL_TOUCHED=TRUE'])
    #os.remove(outRas)
    g = band.ReadAsArray()

    for point in speciesOcc.iterrows():
        xi = int((point[1]['x'] - xMin) / pixelSize)
        yi = int((point[1]['y'] - yMax) / -pixelSize)

        try:
            if g[yi,xi] != noData:
                g[yi,xi] += 1
            else:
                print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
        except:
            print('point ({}, {}) out of bounds'.format(point[1]['x'], point[1]['y']))
            pass


    band.WriteArray(g)
    targetRasSRS = osr.SpatialReference()
    targetRasSRS.ImportFromWkt(srs.ExportToWkt())
    targetRas.SetProjection(targetRasSRS.ExportToWkt())
    band.FlushCache() 
Example 47
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 4 votes vote down vote up
def getValuesAtPoint(indir, rasterfileList, pos, lon, lat, sp):
    #gt(2) and gt(4) coefficients are zero, and the gt(1) is pixel width, and gt(5) is pixel height.
    #The (gt(0),gt(3)) position is the top left corner of the top left pixel of the raster.
    for i, rs in enumerate(rasterfileList):
        print('processing {}'.format(rs))
        presValues = []
        gdata = gdal.Open('{}/{}.tif'.format(indir,rs))
        gt = gdata.GetGeoTransform()
        band = gdata.GetRasterBand(1)
        nodata = band.GetNoDataValue()

        x0, y0 , w , h = gt[0], gt[3], abs(gt[1]), abs(gt[5])
        
        xmin, xmax, ymin, ymax = min(pos[lon]), max(pos[lon]), min(pos[lat]), max(pos[lat])

        # Specify offset and rows and columns to read
        xoff = int((xmin - x0)/w)
        yoff = int((y0 - ymax)/h)
        xcount = int(math.ceil((xmax - xmin)/w)+1)
        ycount = int(math.ceil((ymax - ymin)/h)+1)

        data = band.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float)
        #free memory
        del gdata
        
        if i == 0:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w) - xoff
                Xc = x0 + int((p[1][lon] - x0)/w)*w + w/2 #the cell center x
                y = int((y0 - p[1][lat])/h) - yoff
                Yc = y0 - int((y0 - p[1][lat])/h)*h - h/2 #the cell center y
                try:
                    if data[y,x] != nodata:
                        value = data[y,x]
                    else:
                        value = numpy.nan
                    presVAL = [p[1][sp],p[1][lon],p[1][lat], '{:.6f}'.format(Xc), '{:.6f}'.format(Yc), value]
                    presValues.append(presVAL)
                except:
                    pass
            df = pandas.DataFrame(presValues, columns=['sp', 'x', 'y', 'Xc', 'Yc', rs])
        else:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w) - xoff
                y = int((y0 - p[1][lat])/h) - yoff
                try:
                    if data[y,x] != nodata:
                        presValues.append(data[y,x])
                    else:
                        presValues.append(numpy.nan)
                except:
                    pass
            df[rs] = pandas.Series(presValues)
    del data, band
    print('extracted values written in dataframe')
    return df


#### function to get all pixel center coordinates and corresponding values from rasters 
Example 48
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 4 votes vote down vote up
def getRasterValues(indir, rasterfileList, skipNoData = True):
    
    for i, rs in enumerate(rasterfileList):
        
        if i == 0:
            vList = []
            gdata = gdal.Open('{}/{}.tif'.format(indir,rs))
            gt = gdata.GetGeoTransform()
            band = gdata.GetRasterBand(1)
            nodata = band.GetNoDataValue()
            data = band.ReadAsArray().astype(numpy.float)
            #free memory
            del gdata

            x0, y0 , w , h = gt[0], gt[3], gt[1], gt[5]

            for r, row in enumerate(data):
                x = 0
                for c, column in enumerate(row):
                    if skipNoData == True:
                        if column == nodata:
                            pass
                        else:
                            x = x0 + c*w + w/2
                            y = y0 + r*h + h/2
                            vList.append(['{:.6f}'.format(x),'{:.6f}'.format(y),column])
                    elif skipNoData == False:
                        x = x0 + c*w + w/2
                        y = y0 + r*h + h/2
                        vList.append(['{:.6f}'.format(x),'{:.6f}'.format(y),column])
                        
            df = pandas.DataFrame(vList, columns=['Xc', 'Yc', rs])
            
        else:
            gdata = gdal.Open('{}/{}.tif'.format(indir,rs))
            gt = gdata.GetGeoTransform()
            band = gdata.GetRasterBand(1)
            nodata = band.GetNoDataValue()
            data = band.ReadAsArray().astype(numpy.float)
            #free memory
            del gdata
            if skipNoData == True: 
                vList = [c for r in data for c in r if c != nodata]
            elif skipNoData == False:
                vList = [c for r in data for c in r]
                
            df[rs] = pandas.Series(vList)
    
    del data, band       
    return(df)


# geo raster to numpy array 
Example 49
Project: chorospy   Author: spyrostheodoridis   File: rasterFunc.py    (license) View Source Project 4 votes vote down vote up
def filterByCoverage(vectorFile, rasterFile, covPerc):
    
    srcVector = ogr.Open(vectorFile)
    srcLayer = srcVector.GetLayer()
    # merge all features in one geometry (multi polygone)
    multi  = ogr.Geometry(ogr.wkbMultiPolygon)
    for feature in srcLayer:
        geom = feature.GetGeometryRef()
        multi.AddGeometry(geom)
    
    #attributes of raster file
    rasList = raster2array(rasterFile)

    xsize = rasList[4][0]
    ysize = abs(rasList[4][1])

    pixel_area = xsize*ysize

    rows = rasList[0].shape[0]
    cols = rasList[0].shape[1]

    x1 = rasList[2][0]
    y1 = rasList[2][3]
    
    #iterate over raster cells
    for i in range(rows):
        for j in range(cols):
            ring = ogr.Geometry(ogr.wkbLinearRing)

            ring.AddPoint(x1, y1)
            ring.AddPoint(x1 + xsize, y1)
            ring.AddPoint(x1 + xsize, y1 - ysize)
            ring.AddPoint(x1, y1 - ysize)
            ring.AddPoint(x1, y1)

            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            intersect = multi.Intersection(poly)

            if intersect.ExportToWkt() != 'GEOMETRYCOLLECTION EMPTY':
                perc = (intersect.GetArea()/pixel_area)*100
                if perc > covPerc:
                    rasList[0][i][j] = numpy.nan     
            x1 += xsize
        x1 = rasList[2][0]
        y1 -= ysize
    
    return(rasList[0]) #return the filtered array


# numpy array to geo raster 
Example 50
Project: cv4ag   Author: worldbank   File: ogr2ogr.py    (license) View Source Project 4 votes vote down vote up
def Usage():

    print( "Usage: ogr2ogr [--help-general] [-skipfailures] [-append] [-update] [-gt n]\n" + \
            "               [-select field_list] [-where restricted_where] \n" + \
            "               [-progress] [-sql <sql statement>] \n" + \
            "               [-spat xmin ymin xmax ymax] [-preserve_fid] [-fid FID]\n" + \
            "               [-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]\n" + \
            "               [-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]\n" + \
            "               [-simplify tolerance]\n" + \
            #// "               [-segmentize max_dist] [-fieldTypeToString All|(type1[,type2]*)]\n" + \
            "               [-fieldTypeToString All|(type1[,type2]*)] [-explodecollections] \n" + \
            "               dst_datasource_name src_datasource_name\n" + \
            "               [-lco NAME=VALUE] [-nln name] [-nlt type] [-dim 2|3] [layer [layer ...]]\n" + \
            "\n" + \
            " -f format_name: output file format name, possible values are:")

    for iDriver in range(ogr.GetDriverCount()):
        poDriver = ogr.GetDriver(iDriver)

        if poDriver.TestCapability( ogr.ODrCCreateDataSource ):
            print( "     -f \"" + poDriver.GetName() + "\"" )

    print( " -append: Append to existing layer instead of creating new if it exists\n" + \
            " -overwrite: delete the output layer and recreate it empty\n" + \
            " -update: Open existing output datasource in update mode\n" + \
            " -progress: Display progress on terminal. Only works if input layers have the \"fast feature count\" capability\n" + \
            " -select field_list: Comma-delimited list of fields from input layer to\n" + \
            "                     copy to the new layer (defaults to all)\n" + \
            " -where restricted_where: Attribute query (like SQL WHERE)\n" + \
            " -sql statement: Execute given SQL statement and save result.\n" + \
            " -skipfailures: skip features or layers that fail to convert\n" + \
            " -gt n: group n features per transaction (default 200)\n" + \
            " -spat xmin ymin xmax ymax: spatial query extents\n" + \
            " -simplify tolerance: distance tolerance for simplification.\n" + \
            #//" -segmentize max_dist: maximum distance between 2 nodes.\n" + \
            #//"                       Used to create intermediate points\n" + \
            " -dsco NAME=VALUE: Dataset creation option (format specific)\n" + \
            " -lco  NAME=VALUE: Layer creation option (format specific)\n" + \
            " -nln name: Assign an alternate name to the new layer\n" + \
            " -nlt type: Force a geometry type for new layer.  One of NONE, GEOMETRY,\n" + \
            "      POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION, MULTIPOINT,\n" + \
            "      MULTIPOLYGON, or MULTILINESTRING.  Add \"25D\" for 3D layers.\n" + \
            "      Default is type of source layer.\n" + \
            " -dim dimension: Force the coordinate dimension to the specified value.\n" + \
            " -fieldTypeToString type1,...: Converts fields of specified types to\n" + \
            "      fields of type string in the new layer. Valid types are : \n" + \
            "      Integer, Real, String, Date, Time, DateTime, Binary, IntegerList, RealList,\n" + \
            "      StringList. Special value All can be used to convert all fields to strings.")

    print(" -a_srs srs_def: Assign an output SRS\n"
        " -t_srs srs_def: Reproject/transform to this SRS on output\n"
        " -s_srs srs_def: Override source SRS\n"
        "\n"
        " Srs_def can be a full WKT definition (hard to escape properly),\n"
        " or a well known definition (i.e. EPSG:4326) or a file with a WKT\n"
        " definition." )

    return False