Shapefile creator for modeling domains

Hi all modelers using (or will use) Python!

I start posting some of my Python scripts (likely short ones) that help my routine tasks at my work.

Today’s code is the 1st script. This script generates three shapefiles based on a user-defined domain definition: cell-center, grids, and box.
Output shapefiles will be generated in the same folder that the script is executed. Output files names are predefined based on the input such as domain_name.

Ideally, this code can be improved by adding a function/code lines to convert GRIDDESC files into input formats needed for the script. Either I can add that function if I have time/need or someone else can take my code to do that work.

I hope this is useful for other modelers. This code is not really optimized. If you end up modifying the code for better performance, etc., please share your code with me and other modelers here. Of course, if you find a bug, please let me know in this post too.

================
# -*- coding: utf-8 -*-
"""
This code creates a shapefile using the modeling domain definition given by a user.

by Byeong-Uk Kim
"""
import numpy as np
from collections import OrderedDict
from shapely import speedups
speedups.disable()

from shapely.geometry import mapping, Polygon, Point
import fiona
from fiona.crs import from_string

############################################################################
# INPUT BLOCK
proj4_str = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs"
domain_name = "12US2"
lcc_llx   = -2412000.0 #lower left corner X
lcc_lly   = -1620000.0 #lower left corner Y
cell_size = 12000.0 # in meters
num_cols  = 396
num_rows  = 246
############################################################################


############################################################################
# CODE BLOCK
aqm_lcc_crs = from_string(proj4_str)

# See discussion about the orientation of polygon vertex for the Shapefile
# It is not conformning the OGC standard. :(
#http://gis.stackexchange.com/questions/119150/order-of-polygon-vertices-in-general-gis-clockwise-or-unterclockwise

this_schema = {'geometry': 'Polygon', 'properties': {'COL': 'int', 'ROW': 'int', 'COLROW': 'int'}} # LROW : crrr
with fiona.open("./{0}_grid.shp".format(domain_name), "w", driver="ESRI Shapefile", crs=aqm_lcc_crs, schema=this_schema) as out_shp_file:
    for i_idx in np.arange(num_cols):
        for j_idx in np.arange(num_rows):
            cell_l_x = lcc_llx+i_idx*cell_size
            cell_r_x = lcc_llx+cell_size+i_idx*cell_size
            cell_l_y = lcc_lly+j_idx*cell_size
            cell_u_y = lcc_lly+cell_size+j_idx*cell_size

            poly = Polygon([(cell_l_x, cell_l_y),
                            (cell_l_x, cell_u_y),
                            (cell_r_x, cell_u_y),
                            (cell_r_x, cell_l_y)])
            i = int(i_idx+1)
            j = int(j_idx+1)

            rec = {'geometry': mapping(poly),
                    'id': '-1',
                    'properties': OrderedDict([(u'COL', i), (u'ROW', j), ('COLROW', i*1000+j)]),
                    'type': 'Feature'}
            out_shp_file.write(rec)
out_shp_file.close()

this_schema = {'geometry': 'Point', 'properties': {'COL': 'int', 'ROW': 'int', 'COLROW': 'int'}} # COLROW : crrr
with fiona.open("./{0}_cell_centers.shp".format(domain_name), "w", driver="ESRI Shapefile", s=aqm_lcc_crs, schema=this_schema) as out_shp_file:
    for i_idx in np.arange(num_cols):
        for j_idx in np.arange(num_rows):
            cell_l_x = lcc_llx+i_idx*cell_size
            cell_r_x = lcc_llx+cell_size+i_idx*cell_size
            cell_l_y = lcc_lly+j_idx*cell_size
            cell_u_y = lcc_lly+cell_size+j_idx*cell_size

            cell_x = (cell_l_x+cell_r_x)*0.5
            cell_y = (cell_l_y+cell_u_y)*0.5

            point = Point([cell_x, cell_y])

            i = int(i_idx+1)
            j = int(j_idx+1)

            rec = {'geometry': mapping(point),
                    'id': '-1',
                    'properties': OrderedDict([(u'COL', i), (u'ROW', j), ('COLROW', i*1000+j)]),
                    'type': 'Feature'}
            out_shp_file.write(rec)
out_shp_file.close()

# Adopting dissolve approach from
#http://gis.stackexchange.com/questions/149959/dissolve-polygons-based-on-attributes-with-python-shapely-fiona
from geopandas import GeoSeries, GeoDataFrame

cells = GeoDataFrame.from_file("./{0}_grid.shp".format(domain_name))
cells_geometry = GeoSeries(list(cells.geometry))
cells_dissolved_geometry = cells_geometry.unary_union

# create a geopandas geodataframe, with columns for state and geometry
cells_dissolved = GeoDataFrame(columns=['geometry'], crs=cells.crs)
cells_dissolved["geometry"] =  GeoSeries([cells_dissolved_geometry])
# save to file
cells_dissolved.to_file("./{0}_box.shp".format(domain_name), driver="ESRI Shapefile")
6 Likes

Hi, the script that I generated and use is a little different than yours. It uses the “fauxioapi” wrapper that I posted in this forum to get the grid information from the GRIDDESC file and ogr/osr to do the shapefile. This will create a shapefile with labeled grid cells (column, row, and concatenated col and row as a unique key).

import ogr
import osr
import fauxioapi as io

griddesc = 'GRIDDESC.txt'
grids = ['12US2','36US3']

def draw_grid(grid_name, grid):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    data_source = driver.CreateDataSource('%s.shp' %grid_name)
    srs = osr.SpatialReference()
    srs.ImportFromProj4(grid.proj4())
    with open('%s.pj4' %grid_name, 'w') as pj4:
        pj4.write('%s\n' %grid.proj4())
    layer = data_source.CreateLayer('cells', srs, ogr.wkbPolygon)
    field_name = ogr.FieldDefn('cellid', ogr.OFTString)
    col_len = len(str(grid.NCOLS))
    row_len = len(str(grid.NROWS))
    field_name.SetWidth(col_len+row_len+1)
    layer.CreateField(field_name)
    for col in range(grid.NCOLS):
        for row in range(grid.NROWS):
            feature = ogr.Feature(layer.GetLayerDefn())
            cellid = str(col+1).rjust(col_len,'0') + '!' + str(row+1).rjust(row_len,'0')
            feature.SetField('cellid', cellid)
            poly = []
            x = (col * grid.XCELL) + grid.XORIG
            y = (row * grid.YCELL) + grid.YORIG
            poly.append('%s %s' %(x,y))
            x += grid.XCELL
            poly.append('%s %s' %(x,y))
            y += grid.YCELL
            poly.append('%s %s' %(x,y))
            x -= grid.XCELL
            poly.append('%s %s' %(x,y))
            y -= grid.YCELL
            poly.append('%s %s' %(x,y))
            poly = ogr.CreateGeometryFromWkt('POLYGON((%s))' %','.join(poly))
            feature.SetGeometry(poly)
            layer.CreateFeature(feature)
            feature = None
    data_source = None

for grid_name in grids:
   grid = io.Grid(grid_name, griddesc)
   draw_grid(grid_name, grid)
1 Like

Thank you for sharing your code!

Hi @bukim1,
Thanks for posting this script. I found a typo in the first fiona.open call. hema=this_schema should be schema=this_schema.
There is no typo in the second fiona.open.
Take care,
Sarah

Thank you for identifying the typo. Did that typo cause any error? I am curious.

Yes, it did: fiona.errors.SchemaError: no schema

Thank you. It must be due to copy & paste code lines or something. Maybe it is the best to upload working codes as files to share scripts in the future.

Thank you very much for sharing the code. I am not able to find the fauxioapi. Can you please share it as well?

Hi, I have posted this module to PyPi and github now. You can install it with pip using “pip install fauxioapi”.

Thank you very much for your quick reply. Can I email you for my errors?