search

MapofLife / MOL

Integrating information about species distributions in an effort to support global understanding of the world's biodiversity.
http://mol.org
BSD 3-Clause "New" or "Revised" License
25 stars 2 forks source link

Polygon simplification via PostGIS #91

Closed eightysteele closed 12 years ago

eightysteele commented 12 years ago

Update the workflow README.md with details about how we're simplifying polygons using PostGIS: https://github.com/MapofLife/MOL/blob/develop/workflow/README.md

codyschank commented 12 years ago

I should note Jeremy wrote the shell script for testing ST_SimplifyPerserveTopology, but I can outline it here:

1) import the shapefiles into a postgre datababse using shp2pgsql using flags -D -d -s 4326 http://postgis.refractions.net/docs/ch04.html#shp2pgsql_usage

2) Project the polygons using ST_Transform. We used Web Mercator (SRID 3785) http://www.postgis.org/docs/ST_Transform.html *note: this step is necessary because the ST_SimplifyPreserveTopology only works properly on projected shapes. (from PostGIS In Action, Obe & Hsu 2011) "ST_Simplify and ST_SimplifyPreserveTopology assume planar coordinates. Should youuse these functions with long lat data (SRID 4326), the resultant geometry can range from slightly askew to completely goofy. First transform your long lat to a planar coordinate, apply ST_Simplify, and then transform back to long lat."

3) Simplify the polygons using ST_SimplifyPreserveTopology

4) Export the polygons back to a shapefile (to view in ArcGIS) using pgsql2shp http://postgis.refractions.net/docs/ch04.html#id2627920

eightysteele commented 12 years ago

Good stuff! @codyschank, can you put this into the workflow README. Also, is the shell script that @jmalczyk wrote in the repository?

codyschank commented 12 years ago

I changed the SRID for shp2pgsql to 3857. However, the prj is still not something that ArcGIS can read properly. I wonder how to fix this.

eightysteele commented 12 years ago

What's the shp2pgsql command line you're using?

codyschank commented 12 years ago

shp2pgsql -D -d -s 4326 $file temp | psql -d mol

psql -d mol -c "CREATE TABLE temp_transformed AS SELECT ST_Transform(the_geom, 3857) as the_geom, gid, specid, latin, occcode, date, editsinfo, citation, notes, origin, raster FROM temp"

eightysteele commented 12 years ago

Hmmm. Are you saying that ArcGIS can't read the temp_transformed table? What's the error you're seeing?

codyschank commented 12 years ago

Later on we export back to a shapefile, so ArcGIS is not dealing the database or the temp_transformed table. We did this because we didn't want to have to deal with figuring out how to view a Postgre database in ArcGIS. Anyways, the file opens fine in ArcGIS, but it can't read the prj file, so it doesn't know what the projection is.

This is what the PostGIS processed projection file looks like: PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]

Should look more like this (I think): PROJCS["World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]

jmalczyk commented 12 years ago

This script will be getting an extra line (or two) that reprojects back to WGS84 (EPSG 4326) because that is what CartoDB expects for input data.

On Fri, Feb 17, 2012 at 3:09 PM, codyschank < reply@reply.github.com

wrote:

Later on we export back to a shapefile, so ArcGIS is not dealing the database or the temp_transformed table. We did this because we didn't want to have to deal with figuring out how to view a Postgre database in ArcGIS. Anyways, the file opens fine in ArcGIS, but it can't read the prj file, so it doesn't know what the projection is.

This is what the PostGIS processed projection file looks like: PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]

Should look more like this (I think):

PROJCS["World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]

Reply to this email directly or view it on GitHub: https://github.com/MapofLife/MOL/issues/91#issuecomment-4026985

codyschank commented 12 years ago

Cool. The problem ArcGIS is having reading the prj file is not a critical issue, just a curiosity. I wonder if it would be able to read the prj file if it was in WGS84 though...but we can revisit it later, if we need to.