Database

Notes from CUGOS Holiday Hackfest on 12/21/2012, Allison and Roger

Goal: Simplify the spatial query that compares the input coordinates to the EPSG bounding boxes
Plan: Pre-build bounding box geometry for the EPSG areas and store in database, rather than building them on the fly in the query

Steps:

• Download shapefile with EPSG area boundaries from here: http://www.epsg.org/polygons/
  • Data stored here: /home/projects/projfinder/server/data/EPSG_2012-11-26_shapefile/
• Load EPSG area boundaries into projfinder database:
  • shp2pgsql -s 4326 -W "LATIN1" EPSG_2012-11-26 epsg_poly projfinder > epsg_poly.sql
  • psql -d projfinder -f epsg_poly.sql
• Then create a table with the bounding box:
  • psql -d projfinder
  • create table epsg_poly_bb as SELECT area_code, area_name, ST_Envelope(geom) as geom FROM epsg_poly;
• Update the SQL query used in projfinder.wsgi
  • The bounding box table replaces the epsg_area table
  • The creation of the geometry in the query is no longer needed

New query:

select sp.srid, split_part(sp.srtext,'"',2),  
st_distance(st_transform(st_geometryfromtext('POINT(%s %s)',4326),sp.srid),  
st_geometryfromtext('POINT(%s %s)',sp.srid))   
from spatial_ref_sys as sp,   
epsg_coordinatereferencesystem as cs,   
epsg_poly_bb as bb   
where st_contains(bb.geom, st_geometryfromtext('POINT(%s %s)',4326)) is true   
and cs.area_of_use_code=bb.area_code and   
exists(select 1 from spatial_ref_sys where srid=cs.coord_ref_sys_code)   
and srid=cs.coord_ref_sys_code   
group by sp.srid, sp.auth_name order by st_distance,char_length(split_part(sp.srtext,'"',2)) limit 10;

I tested this using the initial values from the projfinder.com site and compared it with the result from those same values in the original query and they matched.

Here are the two test queries:

select sp.srid, split_part(sp.srtext,'"',2), 
st_distance(st_transform(st_geometryfromtext('POINT(122.1066000000018 48.033000000000044)',4326),sp.srid),
st_geometryfromtext('POINT(1354804 366934)',sp.srid)) 
from spatial_ref_sys as sp, 
epsg_coordinatereferencesystem as cs, 
epsg_poly_bb as bb 
where st_contains(bb.geom, st_geometryfromtext('POINT(-122.1066000000018 48.033000000000044)',4326)) is true 
and cs.area_of_use_code=bb.area_code and exists(select 1 from spatial_ref_sys where srid=cs.coord_ref_sys_code) 
and srid=cs.coord_ref_sys_code and srid=2285 
group by sp.srid, sp.auth_name 
order by st_distance,char_length(split_part(sp.srtext,'"',2)) limit 10;

select sp.srid, split_part(sp.srtext,'"',2), 
st_distance(st_transform(st_geometryfromtext('POINT(122.1066000000018 48.033000000000044)',4326),sp.srid),
st_geometryfromtext('POINT(1354804 366934)',sp.srid)) 
from spatial_ref_sys as sp, epsg_coordinatereferencesystem as cs, 
epsg_area where st_contains(st_geometryfromtext('POLYGON((' || epsg_area.area_west_bound_lon || ' ' || epsg_area.area_south_bound_lat || ',' || epsg_area.area_east_bound_lon || ' ' || epsg_area.area_south_bound_lat || ',' || epsg_area.area_east_bound_lon || ' ' || epsg_area.area_north_bound_lat || ',' || epsg_area.area_west_bound_lon || ' ' || epsg_area.area_north_bound_lat || ',' || epsg_area.area_west_bound_lon || ' ' || epsg_area.area_south_bound_lat || '))',4326),
st_geometryfromtext('POINT(-122.1066000000018 48.033000000000044)',4326)) is true 
and cs.area_of_use_code=epsg_area.area_code and exists(select 1 from spatial_ref_sys where srid=cs.coord_ref_sys_code) 
and srid=cs.coord_ref_sys_code  and srid=2285 
group by sp.srid, sp.auth_name 
order by st_distance,char_length(split_part(sp.srtext,'"',2)) limit 10;

For a more precise comparison, potentially the query could compare to the actual polygon boundaries, rather than the bounding box, but this level of precision may not be needed for this application.