Python Coding of Geospatial Processing in Web-based Mapping Applications

Python has powerful capabilities for coding elements of Web-based mapping applications. This paper highlights examples of analytical geospatial processing services that we have implemented for several Open Source-based development projects, including the Eastern Interconnection States’ Planning Council (EISPC) Energy Zones Mapping Tool (http://eispctools.anl.gov), the Solar Energy Environmental Mapper (http://solarmapper.anl.gov), and the Ecological Risk Calculator (http://bogi.evs.anl.gov/erc/portal). We used common Open Source tools such as GeoServer, PostGIS, GeoExt, and OpenLayers for the basic Web-based portal, then added custom analytical tools to support more advanced functionality. The analytical processes were implemented as Web Processing Services (WPSs) running on PyWPS, a Python implementation of the Open Geospatial Consortium (OGC) WPS. For report tools, areas drawn by the user in the map interface are submitted to a service that utilizes the spatial extensions of PostGIS to generate buffers for use in querying and analyzing the underlying data. Python code then post-processes the results and outputs JavaScript Object Notation (JSON)-formatted data for rendering. We made use of PyWPS’s integration with the Geographic Resources Analysis Support System (GRASS) to implement flexible, user-adjustable suitability models for several renewable energy generation technologies. In this paper, we provide details about the processing methods we used within these project examples.

United States and parts of Canada. The EZMT includes more than 250 map layers, a flexible suitability modeling capability with more than 35 pre-configured models and 65 input modeling layers, and 19 reports that can be run for user-specified areas within the EI. More background about the project is available from [Arg13]. Solar Energy Environmental Mapper (Solar Mapper) [Sol] provides interactive mapping data on utility-scale solar energy resources and related siting factors in the six southwestern states studied in the Solar Energy Development Programmatic Environmental Impact Statement [DOI12]. The application was first launched in December 2010, and a version that has been reengineered with open-source components is scheduled for launch in June 2014. Solar Mapper supports identification and screeninglevel analyses of potential conflicts between development and environmental resources, and is designed primarily for use by regulating agencies, project planners, and public stakeholders. More details about Solar Mapper can be found in [Sol13].
The Ecological Risk Calculator (ERC) [Erc] estimates risk in individual watersheds in the western United States to federally listed threatened and endangered species, and their designated critical habitats from energy-related surface and groundwater withdrawals. The approach takes into account several biogeographical characteristics of watersheds including occupancy, distribution, and imperilment of species, and their sensitivity to impacts from water withdrawals, as well as geophysical characteristics of watersheds known to include designated critical habitats for species of concern. The ERC is intended to help project planners identify potential levels of conflicts related to listed species (and thus the associated regulatory requirements), and is intended to be used as a preliminary screening tool.
Each of these Web-based mapping applications includes both vector (geographic data stored using coordinates) and raster (geographic data stored as a matrix of equally sized cells) spatial data stored in a relational database. For each application, Python was used to add one or more custom geoprocessing, modeling, or reporting services. The following section provides background on the software environment used, followed by specific examples of code with a discussion about the unique details in each.
One of the distinctive elements of geographic data management and processing is the need for coordinate reference systems and coordinate transformations (projections), which are needed to represent areas on the earth's oblate spheroid shape as planar maps and to manage data in Cartesian coordinate systems. These references appear in the code examples as "3857," the European Petroleum Survey Group (EPSG) Spatial Reference ID (SRID) reference for WGS84 Web Mercator (Auxiliary Sphere) and "102003," the USA Contiguous Albers Equal Area Conic projection commonly used for multi-state and national maps of the United States. These standardized EPSG definitions are now maintained by the International Association of Oil & Gas Producers (OGP) Surveying & Positioning Committee [OGP].
The Web Mercator projection has poor properties for many elements of mapping and navigation [NGA] but is used for most current Web-based mapping applications because of the wide availability of high-quality base maps in the Web Mercator projection from providers such as Google Maps. In the Solar Mapper project, we compared area computations in the southwestern United States using Web Mercator against the Albers Equal Area projection and found very large discrepancies in the results (Table  1).
The distortion inherent in world-scale Mercator projections is easily seen by the horizontal expansion of features, which increases dramatically in the higher northern and southern latitudes. In each of our projects, we chose to store local geographic data in Web Mercator to match the base maps and increase performance. However, for geographic processing such as generating buffers and computing lengths and areas, we first convert coordinates to the Albers Equal Area projection to take advantage of the improved properties of that projection.

SOFTWARE ENVIRONMENT
Each of these systems was built with a multi-tier architecture composed of a Javascript/HTML (hypertext markup language) interface built on Bootstrap Many of the examples show geospatial operations using PostGIS. PostGIS extends the functionality of PostgreSQL for raster and vector spatial data with a robust library of functions. We found it to be well documented, reliable, and the "footprint analysis" tools we describe in the examples run significantly faster than similar tools we had previously developed with a popular commerical GIS framework.

EXAMPLES
One of the primary capabilities of each of our Web applications was using an area selected or drawn by the user for analysis (a "footprint"); collecting vector and raster data inside, intersecting, or near the footprint; and compiling it in a report. The first example shows the steps followed through the whole process, including the user interface, and later examples concentrate on refinements of the Python-coded steps.

Full Process for Footprint Analysis of Power Plant Locations Stored as Point Features
This example is from the EZMT and illustrates part of its Power Plant report. The user draws an area of interest over the map ( Figure 1) and specifies other report parameters (Figure 2). The "Launch Report" button submits a request to the Web application  The Web application initiates the report run by making a WPS request to the service, which is implemented in PyWPS. The request is an XML (extensible markup language) document describing the WPS "Execute" operation and is submitted via a hypertext transfer protocol (HTTP) POST. PyWPS receives this POST request, performs some basic validation and preprocessing, and routes the request to the custom WPSProcess implementation for that request. PyWPS then prepares the HTTP response and returns it to the application server. The code below illustrates the major steps used to generate the data for the report.
We use the psycopg2 library to interact with the database, including leveraging the geographic information system (GIS) capabilities of PostGIS.
# Import PostgresSQL library for database queries import psycopg2 The user-specified footprint corresponding to Figure 1 is hardcoded in this example with Web Mercator coordinates specified in meters and using the Well-Known Text (WKT) format.  A database connection is then established, and a cursor is created.
# Make database connection and cursor conn = psycopg2.connect(host=pg_host, database=pg_database, user=pg_user, password=pg_password) cur = self.conn().cursor() Structured Query Language (SQL) is used to (1) convert the Web Mercator footprint to the Albers Equal Area projection, (2) generate a buffer around the Albers version of the footprint, and (3) convert that buffer back to Web Mercator. In these sections, ST_GeomFromText converts WKT to binary geometry, and ST_AsText converts binary geometry back to WKT. Because WKT does not store projection information, it is given as a parameter in ST_GeomFromText.
Once the data have been retrieved, the code compiles it into a Python dictionary which is rendered and returned as a JSON document (excerpt below). This document is retained by the application for eventual rendering into its final form, HTML with the graphs built with ExtJS. Figure 3 shows a portion of the report. "shortname": "power_plant_platts_existing", "feature_type": "point" } }

Footprint Analysis of Transmission Lines Stored as Line Features
Another EISPC report uses a user-specified footprint to analyze electrical transmission line information; however, rather than only listing features inside the footprint as in the previous example, (1) in contrast to points, line features can cross the footprint boundary; and (2) we want to report the total length of the portion within the footprint rather than only listing the matching records. Note that ST_Intersects is used to collect the lines overlapping the footprint, whereas ST_Intersection is used to calculate lengths of only the portion of the lines within the footprint. In addition, the coordinates are transformed into the Albers Equal Area projection for the length computation. sql = "SELECT category, COUNT( * ),sum(ST_Length("+ "ST_Transform(ST_Intersection("+layer+ ".geom,ST_GeomFromText('"+fp_webmerc+ "', 3857)), 102003))) AS sum_length_fp "+ "FROM "+layer+" WHERE ST_Intersects("+layer+ ".geom,ST_GeomFromText('"+fp_webmerc+ "', 3857)) GROUP BY category ORDER BY category" cur.execute(sql) list = [] for row in cur: # Collect results in list of lists...

Footprint Analysis of Watershed Areas Stored as Polygon Features, with Joined Tables
The

Footprint Analysis of Imperiled Species Sensitivity Stored as Raster (Cell-based) Data
Many of the layers used in the mapping tools are stored as raster (cell-based) data rather than vector (coordinate-based) data. The ST_Clip method can retrieve raster or vector data and returns the data within the footprint. The WHERE clause is important for performance because images in the database are usually stored as many records, each with a tile. ST_Intersects restricts the much more processing-intensive ST_Clip method to the tiles overlapping the footprint. When the footprint overlaps multiple image tiles, multiple records are returned to the cursor, and results are combined in the loop. list = [] sql = "SELECT (pvc).value as val,sum((pvc).count) "+ "FROM (SELECT ST_ValueCount(ST_Clip(rast,1, "+ "ST_GeomFromText('"+fp_wkt"', 3857))) as pvc "+ "FROM "+layer+" as x "+ "WHERE ST_Intersects(rast, ST_GeomFromText('"+ fp_wkt"',3857))) as y "+"GROUP BY val ORDER BY val" cur.execute(sql) for row in cur: list "shortname": "imperiled_species_area", "feature_type": "raster" } }

Elevation Profile along User-Specified Corridor Centerline Using Elevation Data Stored as Raster Data
The Corridor Report in the EZMT includes elevation profiles along the user-input corridor centerline. In this example, an elevation layer is sampled along a regular interval along the centerline. First, the coordinate of the sample point is generated with ST_Line_Interpolate_Point, next, the elevation data are retrieved from the layer with ST_Value.

Computation of Suitability for Wind Turbines Using Raster Data Using GRASS
The suitability models implemented in the EZMT use GRASS software for computations, which are accessed in Python through WPSs. The code below shows the main steps followed when running a suitability model in the EZMT. The models use a set of raster layers as inputs, each representing a siting factor such as wind energy level, land cover, environmental sensitivity, proximity to existing transmission infrastructure, etc. Each input layer is coded with values ranging from 0 (Completely unsuitable) to 100 (Completely suitable), and weights are assigned to each layer representing its relative importance. A composite suitability map is computed using a weighted geometric mean. Figure 4 shows the EZMT model launcher with the default settings for land-based wind turbines with 80-meter hub heights. Processing in the Python code follows the same steps that would be used in the command-line interface. First, the processing resolution is set using g.region. Then, the input layers are processed to normalize the weights to sum to 1.0 (this approach simplifies the model computation). Next, an expression is generated, specifying the formula for the model, and r.mapcalc is called to perform the model computation. r.out.gdal is used to export the model result from GRASS format to Geo-Tiff for compatibility with GeoServer, and the projection is set using gdal_translate from the Geospatial Data Abstraction Library [GDAL] plugin for GRASS. # Set the processing resolution WPSProcess.cmd(self, "g.region res=250") outmap = "run"+str(self.  # Run model using r.mapcalc WPSProcess.cmd(self, "r.mapcalc "+outmap+ "="+str(func)) user_dir = "/srv/ez/shared/models/users/"+ str(self.user_id) if not os.path.exists(user_dir): os.makedirs(user_dir) # Export the model result to GeoTIFF format WPSProcess.cmd(self, "r.out.gdal -c input="+ outmap+" output="+outmap+".tif.raw"+ " type=Byte format=GTiff nodata=255 "+ "createopt='TILED=YES', 'BIGTIFF=IF_SAFER'") # Set the projection of the GeoTIFF to EPSG:3857 WPSProcess.cmd(self, "gdal_translate -a_srs EPSG:3857 "+outmap+ ".tif.raw "+user_dir+"/"+outmap+".tif") as ArcGIS and geoDJango, provide frameworks for web-based mapping applications different from the approach we described in the SOFTWARE ENVIRONMENT section. While it is outside the scope of this paper to discuss the merits of these other approaches, we recommend considering them as alternatives when planning projects. The examples in this paper include vector and raster data, as well as code for converting projections, creating buffers, retrieving features within a specified area, computing areas and lengths, computing a raster-based model, and exporting raster results in GeoTIFF format. All examples are written in Python and run within the OGC-compliant WPS framework provided by PyWPS.

CONCLUSIONS
One of the key points we make is that the Web Mercator projection should not be used for generating buffers or computing lengths or areas because of the distortion inherent in the projection. The examples illustrate how these computations can be performed easily in PostGIS. We chose to use the Albers Equal Area projection, which is commonly used for regional and national maps for the United States. Different projections should be used for more localized areas.
So far our Web-based mapping applications include fairly straightforward analysis and modeling services. However, the same approaches can be used for much more sophisticated applications that tap more deeply into PostGIS and GRASS, or the abundant libraries available in the Python ecosystem. Matplotlib, NetworkX, NumPi, RPy2, and SciPy can each be integrated with Python to provide powerful visualization, networking, mathematics, statistical, scientific, and engineering capabilities.

ACKNOWLEDGMENTS
This work was supported by the U.S. Department of Energy, Office of Electricity Delivery and Energy Reliability; and the U.S. Department of Interior, Bureau of Land Management, through U.S. Department of Energy contract DE-AC02-06CH11357. The submitted manuscript has been created by the University of Chicago as Operator of Argonne National Laboratory ("Argonne") under contract No. DE-AC02-06CH11357 with the U.S. Department of Energy. The U.S. Government retains for itself, and others acting on its behalf, a paid-up, nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.