darkblue blog

dmg_83 Okay. I notice that if you a table like CREATE TABLE tbl (x int, y int) and do SELECT (tbl).x, (tbl).y FROM tbl WHERE tbl = ROW(1,2) you’ll indeed get back 1, 2 (assuming such a row exists). That seems a bit weird, since the left side of the expression “tbl = ROW(1,2)” does not seem to be the same as the right, for we can do (tbl).x, but can’t do (ROW(1,2)).x; tbl contains some additional data that ROW(1,2) does not

dmg_83 I’d like to understand the exact rules PG uses for these. I’m guessing it’s similar to something like how python compares a namedtuple and an ordinary tuple (disregarding the names and treating both as tuples for comparison sakes)

RhodiumToad row(1,2) actually uses column names “f1” and “f2” internally but that’s quite difficult to determine (and usually impossible to make use of)
  you can’t do (row(1,2)).f1 for example
  anyway, tbl = row(1,2) works because there’s a generic (record = record) operator defined
  and any composite type can implicitly cast to the “record” pseudotype
  the operator actually looks up the structure of the records at runtime (cached between calls from a single call site)
  the number and types of fields has to be acceptable, but the column names are ignored

RhodiumToad with foo(x) as (values (‘blah blah blah’)) insert into wherever select x from foo;
merlinm or, with vars as (select ‘a’ as a, 2 as b) insert into foo values (select a from vars), select b from vars);
merlinm oops, missing a parent there
merlinm gah
RhodiumToad with vars(a,b) as (values (‘a’, 2)) …
merlinm ah, touche
darkblue_b with .. is a ‘materialized subselect’ under the hood ?
RhodiumToad yes
strk not really materialized, but yes
RhodiumToad it is materialized
strk doesn’t persist,is what I mean
strk not after the statement
darkblue_b ah – its a temporary table… perhaps
RhodiumToad no
darkblue_b ok – well I dont want to waste your time.. just grasping at things.
strk darkblue_b: only exists during the single statement
darkblue_b yes I read that
RhodiumToad it exists only within the scope of the query and it can’t have indexes, but it is materialized in that it is evaluated at most once (lazily)
darkblue_b writes that down

Consider for a moment, a table of records including addresses and data, such that most every address occurs once, but for unknown reasons, also contains more than one duplicate address (656 in 55,279 below). How to remove duplicate addresses? let us count the ways.. Here is one:

The subselect adds one more term, the rank, which is used in the outer expression to write out only the rank one result. Also, notice that SELECT * obviates the need to write all of these column names twice. So this is useful both in the case of true duplication of the entire row, and in other cases where the lines are different in a way that does not pertain to the situation at hand. This example differentiates with an ORDER rank() on one variable to choose which to keep – other methods could be devised.

No doubt someone will say, why dont you do this in one step, in place? The two step here suits me well enough. You can hold the .csv file as intermediary, and also backup the original table if you like.

Compare Two Tables

Given two tables, one of them is somewhat a subset of another, in certain columns..
   e.g.   address,city,state,zip.
How might we determine which entries in table A are not included in table B ? One way is to use EXCEPT.

  select c1,c2,c3,c4 from A except select c1,c2,c3,c4 from B;


An alternative to that is to use WHERE NOT EXISTS. Here is an example using a correlated subquery, and makes use of a rule to shorten the AND comparison sytnax, as follows (I believe both forms run in similar times)

  select c1,c2,c3 from A where not exists
    (select * from B where (B.c1,B.c2,B.c3)=(A.c1,A.c2,A.c3) );


As a first try, before I had seen EXCEPT, I solved this by doing a RIGHT JOIN to the smaller table A, including another field from B and then checking for NULL. This works but was about 10% slower.

thanks to RhodiumToad for guidance

UPDATE via robe2
“I have duplicate city rows in my table , with the same city name and state code but different geopoints, I am storing points for each of the city, I would like to delete cities that have the same data and are within a certain distance from each other, cities with same names but with a reasonable distance I would like to keep.”

  DELETE FROM cities USING cities As ct
    WHERE ct.gid > cities.gid AND ST_DWithin(ct.geom, cities.geom, 0.025) AND cities.name = ct.name;


SYNOPSIS: Using one million randomly generated points, time finding the ten nearest points to the center of the range by using the KNN distance operator versus using a call to PostGIS ST_Distance().

## Use ST_Distance()

explain analyze 
  SELECT ST_Distance(the_geom, ST_GeomFromText('POINT(500 500)')) 
  FROM tgeom ORDER BY ST_Distance(the_geom, ST_GeomFromText('POINT(500 500)')) 
    asc limit 10; 

Total runtime: 1568.499 ms

 QUERY PLAN                                                            
--------------------------------------------------------
 Limit  (cost=289943.64..289943.67 rows=10 width=128) (actual time=1568.452..1568.457 rows=10 loops=1)
   ->  Sort  (cost=289943.64..292443.64 rows=1000000 width=128) (actual time=1568.444..1568.447 rows=10 loops=1)
         Sort Key: (st_distance(the_geom, '01010000000000000000407F400000000000407F40'::geometry))
         Sort Method: top-N heapsort  Memory: 25kB
         ->  Seq Scan on tgeom  (cost=0.00..268334.00 rows=1000000 width=128) (actual time=0.102..1268.725 rows=1000000 loops=1)

## Use KNN Sorting

explain analyze 
  SELECT the_geom < -> ST_GeomFromText('POINT(500 500)') 
  FROM tgeom ORDER BY the_geom < -> ST_GeomFromText('POINT(500 500)') asc limit 10;

Total runtime: 1.847 ms

 QUERY PLAN                                                            
--------------------------------------------------------

 Limit  (cost=0.00..0.74 rows=10 width=128) (actual time=1.476..1.754 rows=10 loops=1)
   ->  Index Scan using tgeom_idx on tgeom  (cost=0.00..74124.46 rows=1000000 width=128) (actual time=1.474..1.748 rows=10 loops=1)
         Order By: (the_geom < -> '01010000000000000000407F400000000000407F40'::geometry)

Complete Setup Instructions for this Example:
Base System: debian 6.0 ‘squeeze’
Add Repositories: backports testing

sudo apt-get -t squeeze-backports install build-essential postgresql-contrib-9.1 postgresql-client-9.1 postgresql-server-dev-9.1 postgresql-9.1 autoconf automake1.9 libtool flex bison gdb

## setup postgres however you are used to working with it..

svn co http://svn.osgeo.org/postgis/branches/2.0/ postgis_20_branch

svn co http://svn.osgeo.org/geos/branches/3.3/ geos_33_branch

svn co http://svn.osgeo.org/metacrs/proj/trunk/proj proj_trunk

sudo apt-get install libexpat1 libexpat1-dev libxslt1.1 libxml2-dev libjson0-dev libcurl4-openssl-dev python2.6-dev

svn checkout http://libkml.googlecode.com/svn/trunk/ libkml_trunk

svn co http://svn.osgeo.org/gdal/branches/1.9/ gdal_19_branch
./configure –with-python –with-curl –without-libtool

## create database testdb ##
psql testdb -f postgis/postgis.sql
psql testdb -f spatial_ref_sys.sql
psql testdb -f raster/rt_pg/rtpostgis.sql

>psql
psql (9.1.6)
Type “help” for help.

testdb=# select postgis_full_version();
POSTGIS=”2.0.2SVN r10643″ GEOS=”3.3.6dev-CAPI-1.7.6″ PROJ=”Rel. 4.8.0, 6 March 2012″ GDAL=”GDAL 1.9.2, released 2012/10/08″ LIBXML=”2.7.8″ LIBJSON=”UNKNOWN” RASTER

CREATE TABLE tgeom (gid serial PRIMARY KEY,the_geom geometry);
INSERT into tgeom(the_geom) select st_geomfromtext('POINT('|| ((random()*1000)::integer)::text || ' ' ||  ((random()*1000)::integer)::text || ')' ) from generate_series(1,1000000);
CREATE INDEX tgeom_idx on tgeom using GIST(the_geom);

## Done Setup

based on
http://www.depesz.com/2010/12/11/waiting-for-9-1-knngist/

UPDATE for debian/ubuntu
   there is now a repository at apt.postgresql.org

SELECT 
  name, state_name, fips, 
  st_summary(the_geom), 
  st_isvalidreason(the_geom) 
FROM usa_counties 
  WHERE not st_isvalid(the_geom);

NOTICE:  Ring Self-intersection at or near point -70.82466 42.26050
-[ RECORD 1 ]----+-------------------------------------------------
name             | Plymouth
state_name       | Massachusetts
fips             | 25023
st_summary       | MultiPolygon[B] with 1 elements
                 |   Polygon[] with 1 rings
                 |    ring 0 has 31 points
st_isvalidreason | Ring Self-intersection[-70.8246609058028 42.2605069326982]

Time: 268.654 ms

UPDATE usa_counties 
  set the_geom = ST_MakeValid(the_geom) where name='Plymouth' and fips = '25023'; 

SELECT 
  name, state_name, fips, 
  st_summary(the_geom), st_isvalidreason(the_geom) 
FROM usa_counties where not st_isvalid(the_geom); 

 name | state_name | fips | st_summary | st_isvalidreason
------+------------+------+------------+------------------ 
(0 rows) 

select postgis_version();
            postgis_version            
---------------------------------------
 2.0 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
(1 row)

  examples of invalid multipolygon

In California, there are multi-county, regional government agencies that perform some functions like planning public transportation systems. The capital city of California is Sacramento, and the regional council is called the Sacramento Area Council of Governments (SACOG).

Using a database of train and light rail stations in a projection where distance measurement is accurate, and the 2010 US Census data, it is possible to find distance from a train station to a census tract or, better yet a census block, like so:

-- return the distance of each transit station from census tract 1768
SELECT 
  b.ogc_fid, 
  st_distance( t.wkb_geometry, b.wkb_geometry)::integer as dist
FROM
  transit_stations_updated t,
  census_2010_tracts b
WHERE
  t.region = 'SACOG' and b.ogc_fid = 1768
ORDER BY
  st_distance( t.wkb_geometry, b.wkb_geometry);

But wait, if we want to get the distance to a transit stop, it is far more useful to find the distance to the nearest transit stop.. and getting distance to one transit stop is straightforward, but what about all transit stops in the district, using one SQL statement for performance?
  For this, we can use the handy WINDOW FUNCTIONS introduced in PostgreSQL 8.4. Let’s break the problem into parts.
  First, list the distances from one census block, in order, using a RANK() function too:
 

-- return the ranked distance of each transit station from census block 433144
SELECT 
  b.ogc_fid as census_block_id,
  t.ogc_fid as transit_station_id,
  st_distance( t.wkb_geometry, b.wkb_geometry)::integer as dist,
  rank() OVER 
   (partition by b.ogc_fid ORDER BY st_distance( t.wkb_geometry, b.wkb_geometry))
FROM
  transit_stations_updated t,
  census_2010_blocks b
WHERE
   t.region = 'SACOG' and b.ogc_fid = 433144
ORDER BY 
   st_distance( t.wkb_geometry, b.wkb_geometry)::integer;

--
 census_block_id | transit_station_id | dist  | rank 
-----------------+--------------------+-------+------
          433144 |                758 |   176 |    1
          433144 |                715 |   369 |    2
          433144 |                716 |  1152 |    3
          433144 |               1137 |  1377 |    4
...

 
Windowing functions have been written to be very general, so the scope of any description is far larger than this post. What we can observe here is that a partition has been declared by census block, and the order within the partition is by distance between census block and transit stop. Now we have almost all of the components in place to solve: for every census block, return only the nearest transit stop. Here I have written a solution using a subselect – there are definitely other ways to write this! Note that a judicious use of indexing can speed things up quite a bit.

-- return only the nearest transit station to each census block and if less than 100 meters
SELECT census_block_id, transit_station_id, dist
FROM
(
  SELECT 
    b.ogc_fid as census_block_id,
    t.ogc_fid as transit_station_id,
    st_distance( t.wkb_geometry, b.wkb_geometry)::integer as dist,
    rank() OVER 
      (partition by b.ogc_fid ORDER BY st_distance( t.wkb_geometry, b.wkb_geometry)) as rank
  FROM
    transit_stations_updated t,
    census_2010_blocks b
  WHERE 
   st_distance(t.wkb_geometry, b.wkb_geometry) < 100 AND
   b.countyfp10 = '067' AND
   t.region = 'SACOG'
) t1

WHERE t1.rank = 1;

Finally, a variation of this problem is to take all locations within a certain buffer distance of at least one transit stop, and record them along with the distance to the nearest transit stop.


----------------------------------------- install ------------
create or replace function r_median(_float8) returns float as '
median(arg1)
' language 'plr';

CREATE AGGREGATE median (
sfunc = plr_array_accum,
basetype = float8,
stype = _float8,
finalfunc = r_median
);



----------------------------------------- test ---------------
select max(a.res),avg(a.res),median(a.res),min(a.res)
from
( select (st_area(wkb_geometry) - shape_area) as res
from landtypes_090911_marin_parcel_adjusted) a;

thuban1 Time: 21255 ms
r900 Time: 13880 ms
ps Time: 9337 ms
macpro Time: 8856 ms
i7-960+ Time: 6980 ms


ps- add standard deviation sd too
pps- just noticed this post