Create grid points

This notebook creates a grid of points across the study area to predict biophony based on the developed model.

import statements



In [1]:

    
from geo.models import Boundary
import numpy
import pandas
import geopandas
import pyproj
from shapely.geometry import Point, Polygon
from geoalchemy2 import Geometry, WKTElement
from sqlalchemy import create_engine

variable declarations



In [2]:

    
database = create_engine('postgresql://user:password@host:port/database')



In [3]:

    
radius = 500



In [4]:

    
spacing = 200



In [5]:

    
width = 15000



In [6]:

    
height = 12000



In [7]:

    
point_crs = pyproj.Proj(init='EPSG:31254')  # MGI / Austria GK West



In [8]:

    
boundary = Boundary.objects.get(name = "study area").geometry.envelope

create grid



In [9]:

    
coords = numpy.array(boundary.coords[0])
bbox = dict()
bbox['left'] = coords[:, 0].min() - ((coords[:, 0].max() - coords[:, 0].min()) / 2)
bbox['bottom'] = coords[:, 1].min() - ((coords[:, 1].max() - coords[:, 1].min()) / 2)
bbox['right'] = bbox['left'] + width
bbox['top'] = bbox['bottom'] + height

rows = int(numpy.ceil((bbox['top'] - bbox['bottom']) / spacing))
columns = int(numpy.ceil((bbox['right'] - bbox['left']) / spacing))

grid = numpy.empty(shape=(rows, columns), dtype=numpy.ndarray)
x_start = bbox['left']
y_start = bbox['bottom']
for row in range(rows):
    for column in range(columns):
        grid[row, column] = [x_start + (spacing * column), y_start + (spacing * row)]

create points



In [10]:

    
# create the points across the grid
points = [Point(p[0], p[1]) for p in grid.ravel()]
# create a dataframe with the point index
ids = pandas.DataFrame({'id':numpy.arange(rows * columns), }).set_index('id')
# create a geodataframe with the index and the points
grid_points = geopandas.GeoDataFrame(ids, crs=point_crs, geometry=points)
# prepare the geometry for the database
grid_points['geom'] = grid_points['geometry'].apply(lambda x: WKTElement(x.wkt, srid=31254))
# drop the geometry column as it is now duplicative
grid_points.drop('geometry', 1, inplace=True)

save grid points to the database



In [11]:

    
grid_points.to_sql(name='geo_grid_points', 
                   con=database, 
                   if_exists='replace', 
                   dtype={'geom':Geometry('POINT', srid=31254)})