Python osgeo.ogr 模块，Open() 实例源码

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

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

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

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

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

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

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

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

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

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()

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)

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()

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)])

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

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()

def _get_ft_ds():
    refresh_token = OAuth2().get_refresh_token()
    ft_driver = ogr.GetDriverByName('GFT')

    return ft_driver.Open('GFT:refresh=' + refresh_token, True)

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

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

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

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

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

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

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

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

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

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

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

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()

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

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 ###

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)

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)

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)

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)

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)

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)

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")

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())

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)

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

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

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

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()

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()

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

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