/*

  SQL statements carrying out geovalidation on vegbien location data.

  Currently hard-coded to operate on the 'geoscrub' table that is

  assumed to have been generated by the geonames scrubbing process

  (see geonames.sql). This table should minimally contain the following

  columns:

     decimallatitude  (latitude, WGS84 decimal degrees)

     decimallongitude (longitude, WGS84 decimal degrees)

     country          (original asserted country name)

     stateprovince    (original asserted stateprovince name)

     county           (original asserted county name)

     countrystd       (GADM2-mapped country name)

     stateprovincestd (GADM2-mapped state/province name)

     countystd        (GADM2-mapped county name)

  (Note that the geonames.sql script also currently creates columns with

  geonames.org IDs for country, stateprovince, and county.)

  After execution of the statements below, the table will contain the

  following new columns (along with geom and geog columns used by

  postgis):

     latlonvalidity        (validity score for lat/lon coordinate)

     countryvalidity       (validity score for country)

     stateprovincevalidity (validity score for stateprovince)

     countyvalidity        (validity score for county)

     pointcountry          (GADM2 country containing the point)

     pointstateprovince    (GADM2 stateprovince containing the point)

     pointcounty           (GADM2 county containing the point)

  Workflow with actual times (as executed 15-Nov-2012)

     12s create point geoms

     16s create point geogs

     23s index geoms

     26s index geogs

     37s index countrystd

     25s index stateprovincestd

      5s vacuum analyze

      1s mark missing coords (UPDATE 0)

    133s mark valid coords (UPDATE 1702989)

      2s mark invalid coords (UPDATE 4981)

     17s vacuum analyze

     27s mark missing countries (UPDATE 222085)

     36s mark non-gadm countries (UPDATE 6547)

   1277s mark point in country (UPDATE 1400627)

   1491s mark point near country (UPDATE 24626)

    573s mark point near country accurate (UPDATE 1182)

      0s mark point near Antarctica (UPDATE 0)

      2s mark point not in country (UPDATE 54085)

      8s vacuum analyze

     15s mark missing stateprovinces (UPDATE 299015)

     12s mark non-gadm stateprovinces (UPDATE 19716)

    920s mark point in stateprovince (UPDATE 1241068)

    864s mark point near stateprovince (UPDATE 35685)

    424s mark point near stateprovince accurate (UPDATE 1051)

      7s mark point not in stateprovince (UPDATE 111282)

       s vacuum analyze

    128s mark missing counties (UPDATE 1431161)

     29s mark non-gadm counties (UPDATE 63014)

     86s mark point in county (UPDATE 174469)

     65s mark point near county (UPDATE 12477)

      ?s mark point near county accurate (UPDATE ?)

      3s mark point not in county (UPDATE 26849)

      ?s vacuum analyze

    890s look up GADM location using given lat/lon (UPDATE 1656709)

  todo:

  * make firm decision about whether we should assume here that *std

    columns contain valid GADM names (or are NULL), or whether we should

    check that here

  Jim Regetz

  NCEAS

  Created Nov 2012

*/

-- generate point geometries from coordinates

SELECT AddGeometryColumn('public', 'geoscrub', 'geom', 4326, 'POINT', 2 );

UPDATE geoscrub

   SET geom = ST_SetSRID(ST_Point(decimallongitude, decimallatitude), 4326);

-- ... create point geographies too ...

-- note that this excludes (leaves null) any occurrences with coordinates

-- that are not valid decimal degrees

ALTER TABLE geoscrub ADD COLUMN geog geography(POINT, 4326);

UPDATE geoscrub

  SET geog = geom::geography

  WHERE geom && ST_MakeBox2D(ST_Point(-180, -90), ST_Point(180, 90));

-- create indexes

CREATE INDEX "geoscrub_geom_gist" ON geoscrub USING GIST (geom);

CREATE INDEX "geoscrub_geog_gist" ON geoscrub USING GIST (geog);

CREATE INDEX "geoscrub_countrystd_idx" ON geoscrub (countrystd);

CREATE INDEX "geoscrub_stateprovincestd_idx" ON geoscrub (stateprovincestd);

CREATE INDEX "geoscrub_countystd_idx" ON geoscrub (countystd);

-- may want to do this first, otherwise the query planner wants to do a

-- hash join instead of nested loop join, and at least when I use a limit

-- statement to time queries on either side of the threshold where it

-- switches, the hash join does waaaaay worse

SET enable_hashjoin = false;

--

-- validate lat/lon coordinates

--

VACUUM ANALYZE geoscrub;

ALTER TABLE geoscrub ADD COLUMN latlonvalidity int;

-- -1 if missing one or both coordinates

UPDATE geoscrub o

    SET latlonvalidity = -1

    WHERE (

        decimallatitude IS NULL

        OR decimallongitude IS NULL

    );

-- 1 if the coordinates falls within the global bounding box

UPDATE geoscrub o

    SET latlonvalidity = 1

    WHERE geom && ST_MakeBox2D(ST_Point(-180, -90), ST_Point(180, 90));

-- 0 otherwise (have coordinates, but not valid)

UPDATE geoscrub o

    SET latlonvalidity = 0

    WHERE o.latlonvalidity IS NULL;

--

-- validate country

--

VACUUM ANALYZE geoscrub;

ALTER TABLE geoscrub ADD COLUMN countryvalidity int;

-- -1 if missing country name

UPDATE geoscrub o

    SET countryvalidity = -1

    WHERE o.country IS NULL;

-- 0 if a country name is asserted, but not recognized among GADM

-- country names

UPDATE geoscrub o

    SET countryvalidity = 0

    WHERE (o.countrystd IS NULL

        OR NOT EXISTS (

           SELECT 1

             FROM gadm2 c

             WHERE o.countrystd=c.name_0

         )

     ) AND o.countryvalidity IS NULL;

-- 3 if the point is in (or on the border of) the asserted country

UPDATE geoscrub AS o

    SET countryvalidity = 3

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND ST_Intersects(o.geom, c.simple_geom)

    ) AND o.countryvalidity IS NULL;

-- 2 for those that aren't 3's, but the point is in within 5km of

-- the asserted country

UPDATE geoscrub o

    SET countryvalidity = 2

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND ST_DWithin(o.geog, c.simple_geog, 5000)

    ) AND o.countryvalidity IS NULL;

-- now recheck the 2's to see if any are within the country based on the

-- original (unsimplified) geometry, and if so, upgrade to 3

UPDATE geoscrub AS o

    SET countryvalidity = 3

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND ST_Intersects(o.geom, c.geom)

    ) AND o.countryvalidity = 2;

-- also handle special case Antartica, which has invalid geography in

-- gadm2 and thus wasn't included in the within-5km check

UPDATE geoscrub AS o

    SET countryvalidity = 3

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND ST_Intersects(o.geom, c.geom)

    ) AND o.countryvalidity < 3

    AND o.countrystd = 'Antarctica';

-- 1 otherwise (have recognized country name, but point is not within

-- 5km of the asserted country)

UPDATE geoscrub o

    SET countryvalidity = 1

    WHERE o.countryvalidity IS NULL;

-----------------------------

-- validate state/province --

-----------------------------

VACUUM ANALYZE geoscrub;

ALTER TABLE geoscrub ADD COLUMN stateprovincevalidity int;

-- -1 if missing stateprovince name

UPDATE geoscrub o

    SET stateprovincevalidity = -1

    WHERE o.country IS NULL

       OR o.stateprovince IS NULL;

-- 0 if a stateprovince name is asserted, but not recognized among GADM

-- stateprovince names

UPDATE geoscrub o

    SET stateprovincevalidity = 0

    WHERE (o.countrystd IS NULL

        OR o.stateprovincestd IS NULL

        OR NOT EXISTS (

           SELECT 1

             FROM gadm2 c

             WHERE o.countrystd=c.name_0

               AND o.stateprovincestd=c.name_1

        )

    ) AND o.stateprovincevalidity IS NULL;

-- 3 if the point is in (or on the border of) the asserted stateprovince

UPDATE geoscrub AS o

    SET stateprovincevalidity = 3

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND o.stateprovincestd=c.name_1

              AND ST_Intersects(o.geom, c.simple_geom)

    ) AND o.countryvalidity = 3

    AND o.stateprovincevalidity IS NULL;

-- for those that aren't 3's, set to 2 if the point is in within 5km of

-- the asserted stateprovince

UPDATE geoscrub o

    SET stateprovincevalidity = 2

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND o.stateprovincestd=c.name_1

              AND ST_DWithin(o.geog, c.simple_geog, 5000)

    ) AND (o.countryvalidity = 2 OR o.countryvalidity = 3)

    AND o.stateprovincevalidity IS NULL;

-- now recheck the 2's to see if any are within the stateprovince based

-- on the original (unsimplified) geometry, and if so, upgrade to 3

UPDATE geoscrub AS o

    SET stateprovincevalidity = 3

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND o.stateprovincestd=c.name_1

              AND ST_Intersects(o.geom, c.geom)

    ) AND o.stateprovincevalidity = 2;

-- 1 otherwise (have recognized stateprovince name, but point is not within

-- 5km of the asserted stateprovince)

UPDATE geoscrub o

    SET stateprovincevalidity = 1

    WHERE o.stateprovincevalidity IS NULL;

-----------------------------

-- validate county/parish --

-----------------------------

ALTER TABLE geoscrub ADD COLUMN countyvalidity int;

-- -1 if missing county name

UPDATE geoscrub o

    SET countyvalidity = -1

    WHERE o.country IS NULL

       OR o.stateprovince IS NULL

       OR o.county IS NULL;

-- 0 if a county name is asserted, but not recognized among GADM

-- county names

UPDATE geoscrub o

    SET countyvalidity = 0

    WHERE (o.countrystd IS NULL

        OR o.stateprovincestd IS NULL

        OR o.countystd IS NULL

        OR NOT EXISTS (

           SELECT 1

             FROM gadm2 c

             WHERE o.countrystd=c.name_0

               AND o.stateprovincestd=c.name_1

               AND o.countystd=c.name_2

        )

    ) AND o.countyvalidity IS NULL;

-- 3 if the point is in (or on the border of) the asserted county

UPDATE geoscrub AS o

    SET countyvalidity = 3

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND o.stateprovincestd=c.name_1

              AND o.countystd=c.name_2

              AND ST_Intersects(o.geom, c.simple_geom)

    ) AND o.stateprovincevalidity = 3

    AND o.countyvalidity IS NULL;

-- for those that aren't 3's, set to 2 if the point is in within 5km of

-- the asserted county

UPDATE geoscrub o

    SET countyvalidity = 2

    WHERE EXISTS (

        SELECT 1

            FROM gadm2 c

            WHERE o.countrystd=c.name_0

              AND o.stateprovincestd=c.name_1

              AND o.countystd=c.name_2

              AND ST_DWithin(o.geog, c.simple_geog, 5000)

    ) AND (o.stateprovincevalidity = 2 OR o.stateprovincevalidity = 3)

    AND o.countyvalidity IS NULL;

-- 1 otherwise (have recognized county name, but point is not within

-- 5km of the asserted county)

UPDATE geoscrub o

    SET countyvalidity = 1

    WHERE o.countyvalidity IS NULL;

--

-- Look up GADM place names based purely on coordinates for all points

--

ALTER TABLE geoscrub ADD COLUMN pointcountry text;

ALTER TABLE geoscrub ADD COLUMN pointstateprovince text;

ALTER TABLE geoscrub ADD COLUMN pointcounty text;

-- note: only using simplified geometry now, for speed reasons

UPDATE geoscrub o

    SET pointcountry = c.name_0,

        pointstateprovince = c.name_1,

        pointcounty = c.name_2

    FROM gadm2 c

    WHERE ST_Intersects(o.geom, c.simple_geom)

      AND o.latlonvalidity = 1;

SET enable_hashjoin = true;