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