davgis.py

# -*- coding: utf-8 -*-
"""
Authors: Gonzalo E. Espinoza-Dávalos
Contact: g.espinoza@un-ihe.org, gespinoza@utexas.edu
Repository: https://github.com/gespinoza/davgis
Module: davgis
Description:
This module is a python wrapper to simplify scripting and automation of common
GIS workflows used in water resources.
"""

from __future__ import division
import os
import math
import tempfile
import warnings
import ogr
import osr
import gdal
import pandas as pd
import netCDF4
from scipy.interpolate import griddata

np = pd.np


def Buffer(input_shp, output_shp, distance):
    """
    Creates a buffer of the input shapefile by a given distance
    """

    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_lyr_defn = inp_lyr.GetLayerDefn()
    inp_srs = inp_lyr.GetSpatialRef()

    # Output
    out_name = os.path.splitext(os.path.basename(output_shp))[0]
    out_driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(output_shp):
        out_driver.DeleteDataSource(output_shp)
    out_source = out_driver.CreateDataSource(output_shp)

    out_lyr = out_source.CreateLayer(out_name, inp_srs, ogr.wkbPolygon)
    out_lyr_defn = out_lyr.GetLayerDefn()

    # Add fields
    for i in range(inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(i)
        out_lyr.CreateField(field_defn)

    # Add features
    for i in range(inp_lyr.GetFeatureCount()):
        feature_inp = inp_lyr.GetNextFeature()
        geometry = feature_inp.geometry()
        feature_out = ogr.Feature(out_lyr_defn)

        for j in range(0, out_lyr_defn.GetFieldCount()):
            feature_out.SetField(out_lyr_defn.GetFieldDefn(j).GetNameRef(),
                                 feature_inp.GetField(j))

        feature_out.SetGeometry(geometry.Buffer(distance))
        out_lyr.CreateFeature(feature_out)
        feature_out = None

    # Save and/or close the data sources
    inp_source = None
    out_source = None

    # Return
    return output_shp


def Feature_to_Raster(input_shp, output_tiff,
                      cellsize, field_name=False, NoData_value=-9999):
    """
    Converts a shapefile into a raster
    """

    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_srs = inp_lyr.GetSpatialRef()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal.GDT_Int16)

    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr,
                            options=["ATTRIBUTE={0}".format(field_name)])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])

    # Save and/or close the data sources
    inp_source = None
    out_source = None

    # Return
    return output_tiff


def List_Fields(input_lyr):
    """
    Lists the field names of input layer
    """
    # Input
    if isinstance(input_lyr, str):
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr, 0)
        inp_lyr = inp_source.GetLayer()
        inp_lyr_defn = inp_lyr.GetLayerDefn()
    elif isinstance(input_lyr, ogr.Layer):
        inp_lyr_defn = input_lyr.GetLayerDefn()

    # List
    names_ls = []

    # Loop
    for j in range(0, inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(j)
        names_ls.append(field_defn.GetName())

    # Save and/or close the data sources
    inp_source = None

    # Return
    return names_ls


def Raster_to_Array(input_tiff, ll_corner, x_ncells, y_ncells,
                    values_type='float32'):
    """
    Loads a raster into a numpy array
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_data_type = inp_band.DataType

    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    ll_x = ll_corner[0]
    ll_y = ll_corner[1]

    top_left_x = ll_x
    top_left_y = ll_y - cellsize_y*y_ncells

    # Change start point
    temp_path = tempfile.mkdtemp()
    temp_driver = gdal.GetDriverByName('GTiff')
    temp_tiff = os.path.join(temp_path, os.path.basename(input_tiff))
    temp_source = temp_driver.Create(temp_tiff, x_ncells, y_ncells,
                                     1, inp_data_type)
    temp_source.GetRasterBand(1).SetNoDataValue(NoData_value)
    temp_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                 top_left_y, rot_2, cellsize_y))
    temp_source.SetProjection(inp_srs)

    # Snap
    gdal.ReprojectImage(inp_lyr, temp_source, inp_srs, inp_srs,
                        gdal.GRA_Bilinear)
    temp_source = None

    # Read array
    d_type = pd.np.dtype(values_type)
    out_lyr = gdal.Open(temp_tiff)
    array = out_lyr.ReadAsArray(0, 0, out_lyr.RasterXSize,
                                out_lyr.RasterYSize).astype(d_type)
    array[pd.np.isclose(array, NoData_value)] = pd.np.nan
    out_lyr = None

    return array


def Resample(input_tiff, output_tiff, cellsize, method=None,
             NoData_value=-9999):
    """
    Resamples a raster to a different spatial resolution
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_data_type = inp_band.DataType

    top_left_x = inp_transform[0]
    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    top_left_y = inp_transform[3]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    # NoData_value = inp_band.GetNoDataValue()

    x_tot_n = inp_lyr.RasterXSize
    y_tot_n = inp_lyr.RasterYSize

    x_ncells = int(math.floor(x_tot_n * (cellsize_x/cellsize)))

    y_ncells = int(math.floor(y_tot_n * (-cellsize_y/cellsize)))

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_source.GetRasterBand(1).SetNoDataValue(NoData_value)
    out_source.SetGeoTransform((top_left_x, cellsize, rot_1,
                                top_left_y, rot_2, -cellsize))
    out_source.SetProjection(inp_srs)

    # Resampling
    method_dict = {'NearestNeighbour': gdal.GRA_NearestNeighbour,
                   'Bilinear': gdal.GRA_Bilinear,
                   'Cubic': gdal.GRA_Cubic,
                   'CubicSpline': gdal.GRA_CubicSpline,
                   'Lanczos': gdal.GRA_Lanczos,
                   'Average': gdal.GRA_Average,
                   'Mode': gdal.GRA_Mode}

    if method in range(6):
        method_sel = method
    elif method in method_dict.keys():
        method_sel = method_dict[method]
    else:
        warnings.warn('Using default interpolation method: Nearest Neighbour')
        method_sel = 0

    gdal.ReprojectImage(inp_lyr, out_source, inp_srs, inp_srs, method_sel)

    # Save and/or close the data sources
    inp_lyr = None
    out_source = None

    # Return
    return output_tiff


def Array_to_Raster(input_array, output_tiff, ll_corner, cellsize,
                    srs_wkt):
    """
    Saves an array into a raster file
    """
    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    y_ncells, x_ncells = input_array.shape
    gdal_datatype = gdaltype_from_dtype(input_array.dtype)

    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal_datatype)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(-9999)

    out_top_left_x = ll_corner[0]
    out_top_left_y = ll_corner[1] + cellsize*y_ncells

    out_source.SetGeoTransform((out_top_left_x, cellsize, 0,
                                out_top_left_y, 0, -cellsize))
    out_source.SetProjection(str(srs_wkt))
    out_band.WriteArray(input_array)

    # Save and/or close the data sources
    out_source = None

    # Return
    return output_tiff


def Clip(input_tiff, output_tiff, bbox):
    """
    Clips a raster given a bounding box
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    top_left_x = inp_transform[0]
    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    top_left_y = inp_transform[3]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    x_tot_n = inp_lyr.RasterXSize
    y_tot_n = inp_lyr.RasterYSize

    # Bounding box
    xmin, ymin, xmax, ymax = bbox

    # Get indices, number of cells, and top left corner
    x1 = max([0, int(math.floor((xmin - top_left_x)/cellsize_x))])
    x2 = min([x_tot_n, int(math.ceil((xmax - top_left_x)/cellsize_x))])
    y1 = max([0, int(math.floor((ymax - top_left_y)/cellsize_y))])
    y2 = min([y_tot_n, int(math.ceil((ymin - top_left_y)/cellsize_y))])

    x_ncells = x2 - x1
    y_ncells = y2 - y1

    out_top_left_x = top_left_x + x1*cellsize_x
    out_top_left_y = top_left_y + y1*cellsize_y

    # Output
    out_array = inp_array[y1:y2, x1:x2]
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform((out_top_left_x, cellsize_x, rot_1,
                                out_top_left_y, rot_2, cellsize_y))
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(out_array)

    # Save and/or close the data sources
    inp_lyr = None
    out_source = None

    # Return
    return output_tiff


def Raster_to_Points(input_tiff, output_shp):
    """
    Converts a raster to a point shapefile
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)

    top_left_x = transform[0]
    cellsize_x = transform[1]
    top_left_y = transform[3]
    cellsize_y = transform[5]
    NoData_value = inp_band.GetNoDataValue()

    x_tot_n = inp_lyr.RasterXSize
    y_tot_n = inp_lyr.RasterYSize

    top_left_x_center = top_left_x + cellsize_x/2.0
    top_left_y_center = top_left_y + cellsize_y/2.0

    # Read array
    array = inp_lyr.ReadAsArray(0, 0, x_tot_n, y_tot_n)  # .astype(pd.np.float)
    array[pd.np.isclose(array, NoData_value)] = pd.np.nan

    # Output
    out_srs = osr.SpatialReference()
    out_srs.ImportFromWkt(inp_srs)
    out_name = os.path.splitext(os.path.basename(output_shp))[0]
    out_driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(output_shp):
        out_driver.DeleteDataSource(output_shp)
    out_source = out_driver.CreateDataSource(output_shp)

    out_lyr = out_source.CreateLayer(out_name, out_srs, ogr.wkbPoint)
    ogr_field_type = ogrtype_from_dtype(array.dtype)
    Add_Field(out_lyr, "RASTERVALU", ogr_field_type)
    out_lyr_defn = out_lyr.GetLayerDefn()

    # Add features
    for xi in range(x_tot_n):
        for yi in range(y_tot_n):
            value = array[yi, xi]
            if ~pd.np.isnan(value):
                feature_out = ogr.Feature(out_lyr_defn)

                feature_out.SetField2(0, value)

                point = ogr.Geometry(ogr.wkbPoint)
                point.AddPoint(top_left_x_center + xi*cellsize_x,
                               top_left_y_center + yi*cellsize_y)

                feature_out.SetGeometry(point)
                out_lyr.CreateFeature(feature_out)

                feature_out = None

    # Save and/or close the data sources
    inp_lyr = None
    out_source = None

    # Return
    return output_shp


def Add_Field(input_lyr, field_name, ogr_field_type):
    """
    Add a field to a layer using the following ogr field types:
    0 = ogr.OFTInteger
    1 = ogr.OFTIntegerList
    2 = ogr.OFTReal
    3 = ogr.OFTRealList
    4 = ogr.OFTString
    5 = ogr.OFTStringList
    6 = ogr.OFTWideString
    7 = ogr.OFTWideStringList
    8 = ogr.OFTBinary
    9 = ogr.OFTDate
    10 = ogr.OFTTime
    11 = ogr.OFTDateTime
    """

    # List fields
    fields_ls = List_Fields(input_lyr)

    # Check if field exist
    if field_name in fields_ls:
        raise Exception('Field: "{0}" already exists'.format(field_name))

    # Create field
    inp_field = ogr.FieldDefn(field_name, ogr_field_type)
    input_lyr.CreateField(inp_field)

    return inp_field


def Spatial_Reference(epsg, return_string=True):
    """
    Obtain a spatial reference from the EPSG parameter
    """
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    if return_string:
        return srs.ExportToWkt()
    else:
        return srs


def List_Datasets(path, ext):
    """
    List the data sets in a folder
    """
    datsets_ls = []
    for f in os.listdir(path):
        if os.path.splitext(f)[1][1:] == ext:
            datsets_ls.append(f)
    return datsets_ls


def NetCDF_to_Raster(input_nc, output_tiff, ras_variable,
                     x_variable='longitude', y_variable='latitude',
                     crs={'variable': 'crs', 'wkt': 'crs_wkt'}, time=None):
    """
    Extract a layer from a netCDF file and save it as a raster file.
    For temporal netcdf files, use the 'time' parameter as:
    t = {'variable': 'time_variable', 'value': '30/06/2017'}
    """
    # Input
    inp_nc = netCDF4.Dataset(input_nc, 'r')
    inp_values = inp_nc.variables[ras_variable]
    x_index = inp_values.dimensions.index(x_variable)
    y_index = inp_values.dimensions.index(y_variable)

    if not time:
        inp_array = inp_values[:]
    else:
        time_variable = time['variable']
        time_value = time['value']
        t_index = inp_values.dimensions.index(time_variable)
        time_index = list(inp_nc.variables[time_variable][:]).index(time_value)
        if t_index == 0:
            inp_array = inp_values[time_index, :, :]
        elif t_index == 1:
            inp_array = inp_values[:, time_index, :]
        elif t_index == 2:
            inp_array = inp_values[:, :, time_index]
        else:
            raise Exception("The array has more dimensions than expected")

    # Transpose array if necessary
    if y_index > x_index:
        inp_array = pd.np.transpose(inp_array)

    # Additional parameters
    gdal_datatype = gdaltype_from_dtype(inp_array.dtype)
    NoData_value = inp_nc.variables[ras_variable]._FillValue

    if type(crs) == str:
        srs_wkt = crs
    else:
        crs_variable = crs['variable']
        crs_wkt = crs['wkt']
        exec('srs_wkt = str(inp_nc.variables["{0}"].{1})'.format(crs_variable,
                                                                 crs_wkt))

    inp_x = inp_nc.variables[x_variable]
    inp_y = inp_nc.variables[y_variable]

    cellsize_x = abs(pd.np.mean([inp_x[i] - inp_x[i-1]
                                 for i in range(1, len(inp_x))]))
    cellsize_y = -abs(pd.np.mean([inp_y[i] - inp_y[i-1]
                                  for i in range(1, len(inp_y))]))

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)

    y_ncells, x_ncells = inp_array.shape

    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal_datatype)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(pd.np.asscalar(NoData_value))

    out_top_left_x = inp_x[0] - cellsize_x/2.0
    if inp_y[-1] > inp_y[0]:
        out_top_left_y = inp_y[-1] - cellsize_y/2.0
        inp_array = pd.np.flipud(inp_array)
    else:
        out_top_left_y = inp_y[0] - cellsize_y/2.0

    out_source.SetGeoTransform((out_top_left_x, cellsize_x, 0,
                                out_top_left_y, 0, cellsize_y))
    out_source.SetProjection(srs_wkt)
    out_band.WriteArray(inp_array)
    out_band.ComputeStatistics(True)

    # Save and/or close the data sources
    inp_nc.close()
    out_source = None

    # Return
    return output_tiff


def Apply_Filter(input_tiff, output_tiff, number_of_passes):
    """
    Smooth a raster by replacing cell value by the average value of the
    surrounding cells
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    top_left_x = inp_transform[0]
    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    top_left_y = inp_transform[3]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Filter
    inp_array[inp_array == NoData_value] = pd.np.nan
    out_array = array_filter(inp_array, number_of_passes)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                top_left_y, rot_2, cellsize_y))
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(out_array)

    # Save and/or close the data sources
    inp_lyr = None
    out_source = None

    # Return
    return output_tiff


def Extract_Band(input_tiff, output_tiff, band_number=1):
    """
    Extract and save a raster band into a new raster
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(band_number)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform(inp_transform)
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(inp_array)

    # Save and/or close the data sources
    inp_lyr = None
    out_source = None

    # Return
    return output_tiff


def Get_Extent(input_lyr):
    """
    Obtain the input layer extent (xmin, ymin, xmax, ymax)
    """
    # Input
    filename, ext = os.path.splitext(input_lyr)
    if ext.lower() == '.shp':
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr)
        inp_lyr = inp_source.GetLayer()
        x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
        inp_lyr = None
        inp_source = None
    elif ext.lower() == '.tif':
        inp_lyr = gdal.Open(input_lyr)
        inp_transform = inp_lyr.GetGeoTransform()
        x_min = inp_transform[0]
        x_max = x_min + inp_transform[1] * inp_lyr.RasterXSize
        y_max = inp_transform[3]
        y_min = y_max + inp_transform[5] * inp_lyr.RasterYSize
        inp_lyr = None
    else:
        raise Exception('The input data type is not recognized')
    return (x_min, y_min, x_max, y_max)


def Interpolation_Default(input_shp, field_name, output_tiff,
                          method='nearest', cellsize=None):
    '''
    Interpolate point data into a raster
    Available methods: 'nearest', 'linear', 'cubic'
    '''
    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_srs = inp_lyr.GetSpatialRef()
    inp_wkt = inp_srs.ExportToWkt()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    ll_corner = [x_min, y_min]
    if not cellsize:
        cellsize = min(x_max - x_min, y_max - y_min)/25.0
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Feature points
    x = []
    y = []
    z = []
    for i in range(inp_lyr.GetFeatureCount()):
        feature_inp = inp_lyr.GetNextFeature()
        point_inp = feature_inp.geometry().GetPoint()

        x.append(point_inp[0])
        y.append(point_inp[1])
        z.append(feature_inp.GetField(field_name))

    x = pd.np.array(x)
    y = pd.np.array(y)
    z = pd.np.array(z)

    # Grid
    X, Y = pd.np.meshgrid(pd.np.linspace(x_min + cellsize/2.0,
                                         x_max - cellsize/2.0,
                                         x_ncells),
                          pd.np.linspace(y_min + cellsize/2.0,
                                         y_max - cellsize/2.0,
                                         y_ncells))

    # Interpolate
    out_array = griddata((x, y), z, (X, Y), method=method)
    out_array = pd.np.flipud(out_array)

    # Save raster
    Array_to_Raster(out_array, output_tiff, ll_corner, cellsize, inp_wkt)

    # Return
    return output_tiff

def get_neighbors(x, y, nx, ny, cells=1):
    """
    Get a list of neighboring cells
    """
    neighbors_ls = [(xi, yi)
                    for xi in range(x - 1 - cells + 1, x + 2 + cells - 1)
                    for yi in range(y - 1 - cells + 1, y + 2 + cells - 1)
                    if (-1 < x <= nx - 1 and -1 < y <= ny - 1 and
                        (x != xi or y != yi) and
                        (0 <= xi <= nx - 1) and (0 <= yi <= ny - 1))]
    return neighbors_ls


def get_mean_neighbors(array, index, include_cell=False):
    """
    Get the mean value of neighboring cells
    """
    xi, yi = index
    nx, ny = array.shape
    stay = True
    cells = 1
    while stay:
        neighbors_ls = get_neighbors(xi, yi, nx, ny, cells)
        if include_cell:
            neighbors_ls = neighbors_ls + [(xi, yi)]
        values_ls = [array[i] for i in neighbors_ls]
        if pd.np.isnan(values_ls).all():
            cells += 1
        else:
            value = pd.np.nanmean(values_ls)
            stay = False
    return value


def array_filter(array, number_of_passes=1):
    """
    Smooth cell values by replacing each cell value by the average value of the
    surrounding cells
    """
    while number_of_passes >= 1:
        ny, nx = array.shape
        arrayf = pd.np.empty(array.shape)
        arrayf[:] = pd.np.nan
        for j in range(ny):
            for i in range(nx):
                arrayf[j, i] = get_mean_neighbors(array, (j, i), True)
        array[:] = arrayf[:]
        number_of_passes -= 1
    return arrayf


def ogrtype_from_dtype(d_type):
    """
    Return the ogr data type from the numpy dtype
    """
    # ogr field type
    if 'float' in d_type.name:
        ogr_data_type = 2
    elif 'int' in d_type.name:
        ogr_data_type = 0
    elif 'string' in d_type.name:
        ogr_data_type = 4
    elif 'bool' in d_type.name:
        ogr_data_type = 8
    else:
        raise Exception('"{0}" is not recognized'.format(d_type))
    return ogr_data_type


def gdaltype_from_dtype(d_type):
    """
    Return the gdal data type from the numpy dtype
    """
    # gdal field type
    if 'int8' == d_type.name:
        gdal_data_type = 1
    elif 'uint16' == d_type.name:
        gdal_data_type = 2
    elif 'int16' == d_type.name:
        gdal_data_type = 3
    elif 'uint32' == d_type.name:
        gdal_data_type = 4
    elif 'int32' == d_type.name:
        gdal_data_type = 5
    elif 'float32' == d_type.name:
        gdal_data_type = 6
    elif 'float64' == d_type.name:
        gdal_data_type = 7
    elif 'bool' in d_type.name:
        gdal_data_type = 1
    elif 'int' in d_type.name:
        gdal_data_type = 5
    elif 'float' in d_type.name:
        gdal_data_type = 7
    elif 'complex' == d_type.name:
        gdal_data_type = 11
    else:
        warnings.warn('"{0}" is not recognized. '
                      '"Unknown" data type used'.format(d_type))
        gdal_data_type = 0
    return gdal_data_type