Re: [postgis-users] postgis-users Digest, Vol 203, Issue 11

2019-01-17 Thread David M. Kaplan

Hi,

I tried adding a gist index to geom::geography and oddly enough it made 
the query a lot slower... Not really sure why, but the best option in my 
case still seems to be to first doing a bounding box reduction in 
lon-lat space and then apply ST_DWithin after casting to geography...


Cheers,
David


On 15/01/2019 21:00, postgis-users-requ...@lists.osgeo.org wrote:
From: Devdatta Tengshe  If you are running a 
query on the Geography, you should create the index on Geography, 
using something like this: Create Index g_geog_idx on mySchema.myTable 
using GIST(geom::geography); Using this increased the speed of my 
query 1000 times. Source: 
https://twitter.com/postgis/status/675001071505383424 Regards, 
Devdatta On Tue, Jan 15, 2019 at 4:07 PM David M. Kaplan 
 wrote:

Hi,

For a while it has not been clear in my head if and when GIST indexes are
used when doing geography based calculations. If the data is stored in a
table with a geometry column that has a GIST index on it, will that index
be used if one does something like
ST_DWithin(a.geom::geography,b.geom,1000)? If not, then if the data was
stored in geography format instead of geometry, would this make the index
more useful?

To give a specific context, I have a table of point geometries with SRID
4326 and a GIST index, and I want to find the number of points in that
table that are within a certain distance of each other point in that table,
but I am not sure what is the most efficient way to do so. To do this, I
have a query of the form:

SELECT a.gid, a.geom, count(*) AS num_points

FROM mytable a

JOIN mytable b

  ON a.gid<>b.gid
  AND ST_DWithin(a.geom::geography,b.geom,1000,FALSE)

GROUP BY a.gid,a.geom

;


This is quite slow, so I presume the GIST index is not being used (and
EXPLAIN didn't show anything that made it clear that it was being used
though bounding box comparisons are included in the join filter). Also
adding a lonlat bounding box comparison does speed the calculation up
significantly:

SELECT a.gid, a.geom, count(*) AS num_points

FROM mytable a

JOIN mytable b

  ON a.gid<>b.gid
  AND a.geom <#> b.geom < 0.02 -- in region where 0.02 degrees is > 1000 m
  AND ST_DWithin(a.geom::geography,b.geom,1000,FALSE)

GROUP BY a.gid,a.geom

;


Is this the best approach or am I missing something? Would using a LATERAL
join improve things?

Thanks,
David


--
**********
David M. Kaplan
Charge de Recherche 1

Institut de Recherche pour le Developpement (IRD)
UMR MARBEC (IRD/Ifremer/CNRS/UMII)
av. Jean Monnet
CS 30171
34203 Sete cedex
France

Email: david.kap...@ird.fr
Phone: +33 (0)4 99 57 32 25
Fax: +33 (0)4 99 57 32 95

http://www.umr-marbec.fr/kaplan-david.html
http://www.davidmkaplan.fr/
**

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] confusion regarding use of indexes when doing geography calculations

2019-01-15 Thread David M. Kaplan

Hi,

For a while it has not been clear in my head if and when GIST indexes 
are used when doing geography based calculations. If the data is stored 
in a table with a geometry column that has a GIST index on it, will that 
index be used if one does something like 
ST_DWithin(a.geom::geography,b.geom,1000)? If not, then if the data was 
stored in geography format instead of geometry, would this make the 
index more useful?


To give a specific context, I have a table of point geometries with SRID 
4326 and a GIST index, and I want to find the number of points in that 
table that are within a certain distance of each other point in that 
table, but I am not sure what is the most efficient way to do so. To do 
this, I have a query of the form:


SELECT a.gid, a.geom, count(*) AS num_points

FROM mytable a

JOIN mytable b

 ON a.gid<>b.gid
 AND ST_DWithin(a.geom::geography,b.geom,1000,FALSE)

GROUP BY a.gid,a.geom

;


This is quite slow, so I presume the GIST index is not being used (and 
EXPLAIN didn't show anything that made it clear that it was being used 
though bounding box comparisons are included in the join filter). Also 
adding a lonlat bounding box comparison does speed the calculation up 
significantly:


SELECT a.gid, a.geom, count(*) AS num_points

FROM mytable a

JOIN mytable b

 ON a.gid<>b.gid
 AND a.geom <#> b.geom < 0.02 -- in region where 0.02 degrees is > 1000 m
 AND ST_DWithin(a.geom::geography,b.geom,1000,FALSE)

GROUP BY a.gid,a.geom

;


Is this the best approach or am I missing something? Would using a 
LATERAL join improve things?


Thanks,
David


--
**********
David M. Kaplan
Charge de Recherche 1

Institut de Recherche pour le Developpement (IRD)
UMR MARBEC (IRD/Ifremer/CNRS/UMII)
av. Jean Monnet
CS 30171
34203 Sete cedex
France

Email: david.kap...@ird.fr
Phone: +33 (0)4 99 57 32 25
Fax: +33 (0)4 99 57 32 95

http://www.umr-marbec.fr/kaplan-david.html
http://www.davidmkaplan.fr/
**

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] very strange bug generating incorrect ST_Distance calculations

2019-01-10 Thread David M. Kaplan
> >/Glad to hear that you seem to have found the issue. If there are 
any other tests I can do to try to diagnose the problem, let me know. /



Can you open an issue in trac (https://trac.osgeo.org/postgis/) with
the details, please? If you don't have an account I can do it myself.



Done! Ticket #4290.




___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] very strange bug generating incorrect ST_Distance calculations

2019-01-10 Thread David M. Kaplan
Glad to hear that you seem to have found the issue. If there are any 
other tests I can do to try to diagnose the problem, let me know.


Cheers,
David

:
You’ve found a bug, it’s in PostGIS (geodetic code is in postgis 
native), and the extra detail that you can get it to manifest with 
MULTIPOLYGON but the same rings in POLYGON don’t manifest is an 
excellent extra detail to help in fixing it.

On 10/01/2019 12:26, David M. Kaplan wrote:

Comparing the results of _ST_DistanceTree and _ST_DistanceUnCached is
a good way to tease out whether there are differences between the
implementations without having to fight the caching machinery along
the way.

P.


Thanks for the responses. I tried with _ST_DistanceTree and 
_ST_DistanceUnCached
and as predicted, tree gives 0 and uncached gives the correct answer. 
What is the appropriate work around to avoid this issue? Who should I 
report this problem to? The maintainers of libgeos?


Also, how should I interpret that this error only occurs when one of 
the geometries is a multipolygon?


In response to Raúl, I found that both queries give 0 for postgis 
2.5.0 r16836, but only the multiline query gives zero for 
POSTGIS="2.4.3 r16312"...


Cheers,
David


On Thu, Jan 10, 2019 at 7:00 AM Raúl Marín Rodríguez
https://lists.osgeo.org/mailman/listinfo/postgis-users>> wrote:
>//>/Hi, />//>/In my case I get dist_geography between gid 1 and gid 3 0 in both 
queries. />//>/This looks like a possible issue around the geos cache, something 
similar />/to https://trac.osgeo.org/postgis/ticket/4269. Can you run them with />/EXPLAIN ANALYZE? />//>/-- />/Raúl Marín Rodríguez />/carto.com />/___ />/postgis-users mailing list />/postgis-users at lists.osgeo.org 
<https://lists.osgeo.org/mailman/listinfo/postgis-users> />/https://lists.osgeo.org/mailman/listinfo/postgis-users /



Hi,

I have found what seems to be an extremely bizarre bug. I am not sure 
if it is a postgis issue or postgresql issue, but it leads 
ST_Distance(geography(Polygon,4326),geography(MultiPolygon,4326)) to 
return 0 (zero) when the true answer is not zero under certain 
peculiar situations.


The table data for a set of polygons that generate the error can be 
downloaded here:


http://www.davidmkaplan.fr/bad_geography_dist_calcs.csv

Once you have that data, the code to generate the error is:

CREATE TEMP TABLE tt (gid int PRIMARY KEY, geom geometry);

\copy tt FROM 'bad_geography_dist_calcs.csv' WITH (FORMAT 
'csv',HEADER TRUE)


-- Incorrect dist_geography for t1.gid=1 and t2.gid=3
SELECT t1.gid AS gid1, t2.gid AS gid2,
   ST_Distance(t1.geom,t2.geom) AS dist_lonlat,
   ST_Distance(ST_Transform(t1.geom,26918),
   ST_Transform(t2.geom,26918)) AS dist_utm,
   ST_Distance(t1.geom::geography,t2.geom::geography) AS 
dist_geography,
ST_Distance((ST_Dump(t1.geom)).geom::geography,(ST_Dump(t2.geom)).geom::geography) 
AS dist_geography_dump

FROM tt t1 JOIN tt t2 ON t1.gid   ST_Distance(t1.geom::geography,t2.geom::geography) AS 
dist_geography,
ST_Distance((ST_Dump(t1.geom)).geom::geography,(ST_Dump(t2.geom)).geom::geography) 
AS dist_geography_dump

FROM tt t1 JOIN tt t2 ON t1.gidHopefully, if I am not crazy, you will have a 0 value of 
dist_geography in the first query for t2.gid=3, but not in the second.


Does anyone have any idea what could be causing this? To be honest, I 
am not really certain that I know where to start do diagnose the 
problem...


The specifics of my installation are:

POSTGIS="2.4.3 r16312" PGSQL="100" GEOS="3.7.0-CAPI-1.11.0 673b9939" 
PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.3.2, released 
2018/09/21" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" RASTER


psql 10.6 (Ubuntu 10.6-0ubuntu0.18.04.1))

Ubuntu 18.04.1 LTS

I have also tested this with postgis 2.5.0 r16836 and there the 
incorrect 0 is in both queries!?!?!?


Thanks,
David




--
**
David M. Kaplan
Charge de Recherche 1

Institut de Recherche pour le Developpement (IRD)
UMR MARBEC (IRD/Ifremer/CNRS/UMII)
av. Jean Monnet
CS 30171
34203 Sete cedex
France

Email:david.kap...@ird.fr
Phone: +33 (0)4 99 57 32 25
Fax: +33 (0)4 99 57 32 95

http://www.umr-marbec.fr/kaplan-david.html
http://www.davidmkaplan.fr/
**


--
**
David M. Kaplan
Charge de Recherche 1

Institut de Recherche pour le Developpement (IRD)
UMR MARBEC (IRD/Ifremer/CNRS/UMII)
av. Jean Monnet
CS 30171
34203 Sete cedex
France

Email: david.kap...@ird.fr
Phone: +33 (0)4 99 57 32 25
Fax: +33 (0)4 99 57 32 95

http://www.umr-marbec.fr/kaplan-david.html
http://www.davidmkaplan.fr/
**

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] very strange bug generating incorrect ST_Distance calculations

2019-01-10 Thread David M. Kaplan


On 10/01/2019 12:26, David M. Kaplan wrote:

Comparing the results of _ST_DistanceTree and _ST_DistanceUnCached is
a good way to tease out whether there are differences between the
implementations without having to fight the caching machinery along
the way.

P.


Thanks for the responses. I tried with _ST_DistanceTree and 
_ST_DistanceUnCached
and as predicted, tree gives 0 and uncached gives the correct answer. 
What is the appropriate work around to avoid this issue? Who should I 
report this problem to? The maintainers of libgeos?


Also, how should I interpret that this error only occurs when one of the 
geometries is a multipolygon?


In response to Raúl, I found that both queries give 0 for postgis 2.5.0 
r16836, but only the multiline query gives zero for POSTGIS="2.4.3 
r16312"...


Cheers,
David



On Thu, Jan 10, 2019 at 7:00 AM Raúl Marín Rodríguez
https://lists.osgeo.org/mailman/listinfo/postgis-users>> wrote:
>//>/Hi, />//>/In my case I get dist_geography between gid 1 and gid 3 0 in both 
queries. />//>/This looks like a possible issue around the geos cache, something 
similar />/to https://trac.osgeo.org/postgis/ticket/4269. Can you run them with />/EXPLAIN ANALYZE? />//>/-- />/Raúl Marín Rodríguez />/carto.com />/___ />/postgis-users mailing list />/postgis-users at lists.osgeo.org 
<https://lists.osgeo.org/mailman/listinfo/postgis-users> />/https://lists.osgeo.org/mailman/listinfo/postgis-users /



Hi,

I have found what seems to be an extremely bizarre bug. I am not sure 
if it is a postgis issue or postgresql issue, but it leads 
ST_Distance(geography(Polygon,4326),geography(MultiPolygon,4326)) to 
return 0 (zero) when the true answer is not zero under certain 
peculiar situations.


The table data for a set of polygons that generate the error can be 
downloaded here:


http://www.davidmkaplan.fr/bad_geography_dist_calcs.csv

Once you have that data, the code to generate the error is:

CREATE TEMP TABLE tt (gid int PRIMARY KEY, geom geometry);

\copy tt FROM 'bad_geography_dist_calcs.csv' WITH (FORMAT 'csv',HEADER 
TRUE)


-- Incorrect dist_geography for t1.gid=1 and t2.gid=3
SELECT t1.gid AS gid1, t2.gid AS gid2,
   ST_Distance(t1.geom,t2.geom) AS dist_lonlat,
   ST_Distance(ST_Transform(t1.geom,26918),
   ST_Transform(t2.geom,26918)) AS dist_utm,
   ST_Distance(t1.geom::geography,t2.geom::geography) AS 
dist_geography,
ST_Distance((ST_Dump(t1.geom)).geom::geography,(ST_Dump(t2.geom)).geom::geography) 
AS dist_geography_dump

FROM tt t1 JOIN tt t2 ON t1.gid   ST_Distance(t1.geom::geography,t2.geom::geography) AS 
dist_geography,
ST_Distance((ST_Dump(t1.geom)).geom::geography,(ST_Dump(t2.geom)).geom::geography) 
AS dist_geography_dump

FROM tt t1 JOIN tt t2 ON t1.gidHopefully, if I am not crazy, you will have a 0 value of 
dist_geography in the first query for t2.gid=3, but not in the second.


Does anyone have any idea what could be causing this? To be honest, I 
am not really certain that I know where to start do diagnose the 
problem...


The specifics of my installation are:

POSTGIS="2.4.3 r16312" PGSQL="100" GEOS="3.7.0-CAPI-1.11.0 673b9939" 
PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.3.2, released 
2018/09/21" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" RASTER


psql 10.6 (Ubuntu 10.6-0ubuntu0.18.04.1))

Ubuntu 18.04.1 LTS

I have also tested this with postgis 2.5.0 r16836 and there the 
incorrect 0 is in both queries!?!?!?


Thanks,
David




--
**
David M. Kaplan
Charge de Recherche 1

Institut de Recherche pour le Developpement (IRD)
UMR MARBEC (IRD/Ifremer/CNRS/UMII)
av. Jean Monnet
CS 30171
34203 Sete cedex
France

Email: david.kap...@ird.fr
Phone: +33 (0)4 99 57 32 25
Fax: +33 (0)4 99 57 32 95

http://www.umr-marbec.fr/kaplan-david.html
http://www.davidmkaplan.fr/
**

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] very strange bug generating incorrect ST_Distance calculations

2019-01-10 Thread David M. Kaplan

Hi,

I have found what seems to be an extremely bizarre bug. I am not sure if 
it is a postgis issue or postgresql issue, but it leads 
ST_Distance(geography(Polygon,4326),geography(MultiPolygon,4326)) to 
return 0 (zero) when the true answer is not zero under certain peculiar 
situations.


The table data for a set of polygons that generate the error can be 
downloaded here:


http://www.davidmkaplan.fr/bad_geography_dist_calcs.csv

Once you have that data, the code to generate the error is:

CREATE TEMP TABLE tt (gid int PRIMARY KEY, geom geometry);

\copy tt FROM 'bad_geography_dist_calcs.csv' WITH (FORMAT 'csv',HEADER TRUE)

-- Incorrect dist_geography for t1.gid=1 and t2.gid=3
SELECT t1.gid AS gid1, t2.gid AS gid2,
   ST_Distance(t1.geom,t2.geom) AS dist_lonlat,
   ST_Distance(ST_Transform(t1.geom,26918),
   ST_Transform(t2.geom,26918)) AS dist_utm,
   ST_Distance(t1.geom::geography,t2.geom::geography) AS 
dist_geography,
ST_Distance((ST_Dump(t1.geom)).geom::geography,(ST_Dump(t2.geom)).geom::geography) 
AS dist_geography_dump

FROM tt t1 JOIN tt t2 ON t1.gid   ST_Distance(t1.geom::geography,t2.geom::geography) AS 
dist_geography,
ST_Distance((ST_Dump(t1.geom)).geom::geography,(ST_Dump(t2.geom)).geom::geography) 
AS dist_geography_dump

FROM tt t1 JOIN tt t2 ON t1.gidHopefully, if I am not crazy, you will have a 0 value of dist_geography 
in the first query for t2.gid=3, but not in the second.


Does anyone have any idea what could be causing this? To be honest, I am 
not really certain that I know where to start do diagnose the problem...


The specifics of my installation are:

POSTGIS="2.4.3 r16312" PGSQL="100" GEOS="3.7.0-CAPI-1.11.0 673b9939" 
PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.3.2, released 2018/09/21" 
LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.2.1" RASTER


psql 10.6 (Ubuntu 10.6-0ubuntu0.18.04.1))

Ubuntu 18.04.1 LTS

I have also tested this with postgis 2.5.0 r16836 and there the 
incorrect 0 is in both queries!?!?!?


Thanks,
David


--
**
David M. Kaplan
Charge de Recherche 1

Institut de Recherche pour le Developpement (IRD)
UMR MARBEC (IRD/Ifremer/CNRS/UMII)
av. Jean Monnet
CS 30171
34203 Sete cedex
France

Email: david.kap...@ird.fr
Phone: +33 (0)4 99 57 32 25
Fax: +33 (0)4 99 57 32 95

http://www.umr-marbec.fr/kaplan-david.html
http://www.davidmkaplan.fr/
**

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Aggregating rasters by adding and other confusions

2018-03-12 Thread David M. Kaplan
n place.
> 
> What it seems to do is two passes
> 
>  
> 
> 1 pass does an ST_Union using 'COUNT'  (so in your case you'd get
> numbers between 0 and 3, 0 being no rasters considered having any
> data)
> 
> 2nd pass does an ST_Union using 'SUM'  
> 
> And returns a new raster  where each pixel is  SUM/COUNT
> 
>  
> 
> Basic code is here:
> 
>  
> 
> http://postgis.net/docs/doxygen/2.4/da/dde/rtpg__mapalgebra_8c_a1d940
> 65e6cef5d5d61417b82b2cf4fb6.html#a1d94065e6cef5d5d61417b82b2cf4fb6
> 
> (note it only computes a value if the first band (which I presume to
> be the count band)  value > 0  and has no nodata value. )
> 
> I don't know why a count band would have a nodata value aside from
> when it's 0
> 
>  
> 
> What does COUNT return for you?  
> 
>  
> 
> Thanks,
> 
> Regina
> 
>  
> 
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of David M. Kaplan
> Sent: Tuesday, March 06, 2018 10:29 AM
> To: postgis-users@lists.osgeo.org
> Subject: Re: [postgis-users] Aggregating rasters by adding and other
> confusions
> 
>  
> 
> Hi,
> 
>  
> 
> Yes, this does seem to be the solution. For completeness, the
> following did the same thing as my aggregate function:
> 
>  
> 
> SELECT ST_Union(rast,1,'SUM') FROM blah;
> 
> This includes replacing 0's with NULL's - ST_Union had the same
> behavior at my aggregate function. 
> 
>  
> 
> Can you explain this? For my sum, this isn't a big deal, but
> ST_Union(... 'MEAN') seems to be treating 0's at NULL values, thereby
> removing them from the mean calculation:
> 
>  
> 
> db=# SELECT (ST_SummaryStats(ST_Union(rast,1,'SUM'))).sum FROM blah;
> sum  
> --
>  2792
> (1 row)
>  
> db=# SELECT (ST_SummaryStats(ST_Union(rast,1,'MEAN'))).sum FROM blah;
>  sum  
> --
>  1346
> (1 row)
>  
> 
> The table "blah" has 3 rasters with no null values, so I would expect
> the sum of the mean to be 1/3 the sum of the sum, but clearly this is
> not the case.
> 
>  
> 
> I tried playing around with the no data value for my rasters to see
> if this would change things. It does change things, but not in any
> way I can explain:
> 
>  
> 
> fads=# SELECT
> (ST_SummaryStats(ST_Union(ST_SetBandNoDataValue(rast,1,0),1,'MEAN')))
> .sum FROM blah;
>  sum  
> --
>  2025
> (1 row)
>  
> fads=# SELECT
> (ST_SummaryStats(ST_Union(ST_SetBandNoDataValue(rast,1,-
> 999),1,'MEAN'))).sum FROM blah;
>  sum  
> --
>  2025
> (1 row)
>  
> 
> Setting the no data value has no effect on the sum of the sum.
> 
>  
> 
> Does any of this make sense?
> 
>  
> 
> Thanks,
> 
> David
> 
>  
> 
>  
> 
> On Mon, 2018-03-05 at 12:00 -0800, postgis-users-request@lists.osgeo.
> org <mailto:postgis-users-requ...@lists.osgeo.org>  wrote:
> 
> Date: Mon, 5 Mar 2018 12:05:40 -0500
> From: "Regina Obe" <l...@pcorp.us <mailto:l...@pcorp.us> >
> To: "'PostGIS Users Discussion'" <postgis-users@lists.osgeo.org  to:postgis-users@lists.osgeo.org> >
> Subject: Re: [postgis-users] Aggregating rasters by adding and other
> confusions
> Message-ID: <001f01d3b4a4$30d1fa20$9275ee60$@pcorp.us <mailto:001f01d
> 3b4a4$30d1fa20$9275ee60$@pcorp.us> >
> Content-Type: text/plain; charset="utf-8"
>  
> I think what you are looking for is ST_Union (.. SUM)  note this has
> union types – FIRST, MIN, MAX, COUNT, SUM, MEAN, RANGE
>  
>  
>  
> http://postgis.net/docs/manual-2.4/RT_ST_Union.html
>  
>  
>  
>  
>  
>  
>  
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of David M. Kaplan
> Sent: Monday, March 05, 2018 10:03 AM
> To: postgis-users@lists.osgeo.org <mailto:postgis-users@lists.osgeo.o
> rg> 
> Subject: [postgis-users] Aggregating rasters by adding and other
> confusions
>  
>  
>  
> Hi,
>  
>  
>  
> I have recently started working with the postgis raster
> functionality. In general, I have found this really useful and have
> been able to do some neat things fairly simply with this raster
> functionality. Nevertheless, there are a few basic things that I am
> confused about and I was hoping someone could give me a hand.
>  
>  
>  
> (1) First of all, I have a table with a bunch of rasters that have
> the same extent, alignment, scale, etc. and I want to aggregate them
> together into a single raster using pixel-by-pixel addition. It seems
> like there should be a function to do this, but I can't find one. Is
> there an aggregate

Re: [postgis-users] Aggregating rasters by adding and other confusions

2018-03-06 Thread David M. Kaplan
Hi,

Yes, this does seem to be the solution. For completeness, the following
did the same thing as my aggregate function:

SELECT ST_Union(rast,1,'SUM') FROM blah;
 
This includes replacing 0's with NULL's - ST_Union had the same
behavior at my aggregate function. 

Can you explain this? For my sum, this isn't a big deal, but
ST_Union(... 'MEAN') seems to be treating 0's at NULL values, thereby
removing them from the mean calculation:

db=# SELECT (ST_SummaryStats(ST_Union(rast,1,'SUM'))).sum FROM blah; sum  
--
 2792
(1 row)

db=# SELECT (ST_SummaryStats(ST_Union(rast,1,'MEAN'))).sum FROM blah;
 sum  
--
 1346
(1 row)

The table "blah" has 3 rasters with no null values, so I would expect
the sum of the mean to be 1/3 the sum of the sum, but clearly this is
not the case.

I tried playing around with the no data value for my rasters to see if
this would change things. It does change things, but not in any way I
can explain:

fads=# SELECT 
(ST_SummaryStats(ST_Union(ST_SetBandNoDataValue(rast,1,0),1,'MEAN'))).sum FROM 
blah;
 sum  
--
 2025
(1 row)

fads=# SELECT 
(ST_SummaryStats(ST_Union(ST_SetBandNoDataValue(rast,1,-999),1,'MEAN'))).sum 
FROM blah;
 sum  
--
 2025
(1 row)

Setting the no data value has no effect on the sum of the sum.

Does any of this make sense?

Thanks,
David


On Mon, 2018-03-05 at 12:00 -0800, postgis-users-
requ...@lists.osgeo.org wrote:
> Date: Mon, 5 Mar 2018 12:05:40 -0500
> From: "Regina Obe" <l...@pcorp.us>
> To: "'PostGIS Users Discussion'" <postgis-users@lists.osgeo.org>
> Subject: Re: [postgis-users] Aggregating rasters by adding and other
> confusions
> Message-ID: <001f01d3b4a4$30d1fa20$9275ee60$@pcorp.us>
> Content-Type: text/plain; charset="utf-8"
> 
> I think what you are looking for is ST_Union (.. SUM)  note this has union 
> types – FIRST, MIN, MAX, COUNT, SUM, MEAN, RANGE
> 
>  
> 
> http://postgis.net/docs/manual-2.4/RT_ST_Union.html
> 
>  
> 
>  
> 
>  
> 
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf 
> Of David M. Kaplan
> Sent: Monday, March 05, 2018 10:03 AM
> To: postgis-users@lists.osgeo.org
> Subject: [postgis-users] Aggregating rasters by adding and other confusions
> 
>  
> 
> Hi,
> 
>  
> 
> I have recently started working with the postgis raster functionality. In 
> general, I have found this really useful and have been able to do some neat 
> things fairly simply with this raster functionality. Nevertheless, there are 
> a few basic things that I am confused about and I was hoping someone could 
> give me a hand.
> 
>  
> 
> (1) First of all, I have a table with a bunch of rasters that have the same 
> extent, alignment, scale, etc. and I want to aggregate them together into a 
> single raster using pixel-by-pixel addition. It seems like there should be a 
> function to do this, but I can't find one. Is there an aggregate 
> "ST_MapAlgebra" function? 
> 
>  
> 
> Given that I couldn't find one, I defined an aggregate function as follows:
> 
>  
> 
> CREATE OR REPLACE FUNCTION AddRasters(r1 raster, r2 raster)
>RETURNS raster AS
> $BODY$
> SELECT ST_MapAlgebra($1,$2,'[rast1]+[rast2]');
> $BODY$ LANGUAGE 'sql' IMMUTABLE STRICT;
>  
> CREATE AGGREGATE sum (raster)
> (
> sfunc = AddRasters,
> stype = raster
> );
>  
> 
> (2) This seems to work, but it has the unexpected behavior that it replaces 0 
> values with NULL. In my case, this is fine, but I am wondering why it does 
> this? I can't find anything that indicates that it should be replacing zeros 
> with NULL. Here is the metadata associated with one of my rasters (the others 
> are similar):
> 
>  
> 
> # SELECT ST_BandMetadata(rast), ST_Metadata(rast), ST_SummaryStats(rast) FROM 
> blah;
> -[ RECORD 1 ]---+---
> st_bandmetadata | (16BUI,,f,)
> st_metadata | (-180,90,360,180,1,-1,0,0,4326,1)
> st_summarystats | (64800,417,0.00643518518518519,0.223617719977485,0,46)
>  
> 
> I have not defined a nodataval for these layers and the original layers have 
> no NULL values.
> 
>  
> 
> (3) Is there a postgis command to turn all the NULL values back into zeros?
> 
>  
> 
> (4) I was also considering just defining the '+' operator for raster + raster 
> to be pixel-by-pixel addition. Is there any reason that I wouldn't want to do 
> this?
> 
>  
> 
> (5) Finally, I have been visualizing results with QGIS using the DB Manager. 
> However, I don't see how to select a row from a raster table and incorporate 
> just that row into the canvas. Is there a way to do this?
> 
>  
> 
> Thanks for the assistance,
> 
> David___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Aggregating rasters by adding and other confusions

2018-03-05 Thread David M. Kaplan
Hi,

I have recently started working with the postgis raster functionality.
In general, I have found this really useful and have been able to do
some neat things fairly simply with this raster functionality.
Nevertheless, there are a few basic things that I am confused about and
I was hoping someone could give me a hand.

(1) First of all, I have a table with a bunch of rasters that have the
same extent, alignment, scale, etc. and I want to aggregate them
together into a single raster using pixel-by-pixel addition. It seems
like there should be a function to do this, but I can't find one. Is
there an aggregate "ST_MapAlgebra" function? 

Given that I couldn't find one, I defined an aggregate function as
follows:

CREATE OR REPLACE FUNCTION AddRasters(r1 raster, r2
raster)   RETURNS raster AS$BODY$SELECT
ST_MapAlgebra($1,$2,'[rast1]+[rast2]');$BODY$ LANGUAGE 'sql' IMMUTABLE
STRICT;
CREATE AGGREGATE sum (raster)(sfunc = AddRasters,stype =
raster);
(2) This seems to work, but it has the unexpected behavior that it
replaces 0 values with NULL. In my case, this is fine, but I am
wondering why it does this? I can't find anything that indicates that
it should be replacing zeros with NULL. Here is the metadata associated
with one of my rasters (the others are similar):

# SELECT ST_BandMetadata(rast), ST_Metadata(rast),
ST_SummaryStats(rast) FROM blah;-[ RECORD 1 ]---+
---st_bandmetadata |
(16BUI,,f,)st_metadata | (-180,90,360,180,1,-
1,0,0,4326,1)st_summarystats |
(64800,417,0.00643518518518519,0.223617719977485,0,46)
I have not defined a nodataval for these layers and the original layers
have no NULL values.

(3) Is there a postgis command to turn all the NULL values back into
zeros?

(4) I was also considering just defining the '+' operator for raster +
raster to be pixel-by-pixel addition. Is there any reason that I
wouldn't want to do this?

(5) Finally, I have been visualizing results with QGIS using the DB
Manager. However, I don't see how to select a row from a raster table
and incorporate just that row into the canvas. Is there a way to do
this?

Thanks for the assistance,
David




-- 
**David M. KaplanCharge de Recherche 1
Institut de Recherche pour le Developpement (IRD)UMR MARBEC 
(IRD/Ifremer/CNRS/UMII)av. Jean MonnetCS 3017134203 Sete cedexFrance
Email: david.kaplan@ird.frPhone: +33 (0)4 99 57 32 25Fax: +33 (0)4 99 57 32 95
http://www.umr-marbec.fr/kaplan-david.htmlhttp://www.davidmkaplan.fr/**___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] error while loading dat for tiger geocoder

2017-12-19 Thread David M. Kaplan

Hi,

I am in the process of setting up the tiger geocoder for use with 
addresses from the state of virginia. I am using PostgreSQL 9.3.20, 
postgis 2.2.2 and postgis_tiger_geocoder 2.2.2.


The state loader script (for linux shell) generates an SQL error towards 
the end of the script relating to the following line (line 210 in my 
script):


${PSQL} -c "ALTER TABLE tiger_staging.VA_tabblock RENAME geoid10 TO 
tabblock_id;"


The error it generates is:

ERROR:  relation "tiger_staging.va_tabblock" does not exist

I am under the impression that this error is harmless as at this stage 
the tiger_staging schema is empty because the function 
loader_load_staged_data() has already been run, but I want to make sure.


If this is indeed a harmless error, then why is this line in the script?

Thanks,
David

--
**********
David M. Kaplan
Charge de Recherche 1

Institut de Recherche pour le Developpement (IRD)
UMR MARBEC (IRD/Ifremer/CNRS/UMII)
av. Jean Monnet
CS 30171
34203 Sete cedex
France

Phone: +33 (0)4 99 57 32 25
Fax: +33 (0)4 99 57 32 95

http://www.umr-marbec.fr/kaplan-david.html
**

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users