Re: [postgis-users] PostGIS verse those in Azure Cloud

[David Haynes]

Suprio Ray and Ahmed Eldway are two researchers that do a lot of spatial
computation work on big data (Hadoop, Spark, and parallel computation
frameworks). In my opinion, PostGIS is a more robust tool for spatial
operations. However other tools and platforms can be very good for specific
spatial operations.

On Tue, Aug 29, 2023, 12:40 PM Shaozhong SHI  wrote:

> Geospatial capability has varied maturity in different systems, particular
> of interest is those in Azure Cloud as compared to PostGIS.
>
> Is there any publication on this topic?
>
> Regards,
>
> David
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Raster database design

[David Haynes]

It would make more sense to use a view than another table. If someone
forgets the where f.raster_name check the query time will increase
dramatically.

On Sun, Jul 2, 2023 at 7:40 PM Simon SPDBA Greener 
wrote:

> Currently when I load a raster (GDAL) into PostGIS I get a table per
> raster and a row per tile (if tiled).
>
> Is it possible to store multiple rasters within the same table, with a
> descriminator column to identify which rows belong to which raster?
>
> Such storage makes queries like this possible.
>
> SELECT ST_SummaryStats(rast)
>FROM flood as f
> WHERE f.raster_name = 'toddriver_2019_3p0m';
>
> Or perhaps a better question: what are the pros/cons of doing so?
>
> create table toddriver_2019_2p5m (
> rid integer,
> raster_name text,
> rast raster );
>
> I would assume that one could not apply the constraints to this table.
>
> regards
>
> Simon
>
> --
> Simon Greener
> 39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia
> (m) +61 418 396 391
> (w) www.spdba.com.au
> (m) si...@spdba.com.au
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] parallel-processing multiple similar query tasks - any example?

[David Haynes]

I think you need to provide more information as parallel queries in a
database can be done in many ways. You could look at CitusDB, which
supports PostGIS if you are trying to do parallel queries. PostgreSQL also
supports workers.
https://www.postgresql.org/docs/current/how-parallel-query-works.html

On Wed, Apr 27, 2022 at 6:33 PM Shaozhong SHI 
wrote:

> multiple similar query tasks are as follows:
>
> select * from a_table where country ='UK'
> select * from a_table where country='France'
> and so on
>
> How best to parallel-processing such types of multiple similar query tasks?
>
> Any example available?
>
> Regards,
>
> David
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Enabling output GDAL Driver

[David Haynes]

Hello List,

Running into an error running a 3rd party application GeoServer (2.16) to
visualize rasters stored in PostgreSQL. It seems my PostgreSQL DB has not
enabled the output GDAL driver.

ERROR: rt_raster_to_gdal: Could not load the output GDAL driver

2022-03-26 10:32:24,110 ERROR [jdbc.custom] - ERROR: rt_raster_to_gdal:
Could not load the output GDAL driver
  Where: PL/pgSQL function st_astiff(raster,text[],integer) line 24 at
RETURN
org.postgresql.util.PSQLException: ERROR: rt_raster_to_gdal: Could not load
the output GDAL driver

Overstack suggested this, but it hasn't worked. Do I have to recompile?

SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Geometry Viewer in PgAdmin

[David Haynes]

Here is a github repository from my class that walks you through how to
connect PostgreSQL to QGIS
https://github.com/MGIS-UMN/GIS5577_week5

On Wed, Mar 9, 2022 at 12:43 PM Bo Victor Thomsen <
bo.victor.thom...@gmail.com> wrote:

> Actually, you can upload basically any QGIS vector based layer directly
> into a PostGIS database (connected to QGIS). No need for a shapefile detour
>
> Med venlig hilsen / Kind regards
>
> Bo Victor Thomsen
>
> Den 09-03-2022 kl. 18:45 skrev ruve...@beamerbrooks.com:
>
> Yet another reason to install QGIS is that it will read just about any GIS
> data format; the layer can then be saved as a shapefile which can be
> imported into PostGIS using the PostGIS Shapefile import tool.
>
> Ruven Brooks
>
> On 3/8/2022 9:59 PM, Brent Wood wrote:
>
> I'd also look as DBeaver, it has a free Community Edition and a nice (but
> basic) geometry viewer for Postgis. Built in support for a good set of
> backgrounds as well. Postgis columns are auto displayed as WKT which I find
> useful.
>
> For a user (vs admin) of a Postgis DB I think DBeaver is nicer than
> PgAdmin. It does primitive ERD's, which I don't think PgAdmin can do yet,
> but that might have changed since I last looked...
>
> And you can access many other databases than just Postgres.
>
> But I agree with Ruven, if you want to do more that a casual perusal of
> Postgis data, go QGIS.
>
>
> Brent Wood
>
>
> On Wednesday, March 9, 2022, 03:41:01 PM GMT+13, ruve...@beamerbrooks.com
>   wrote:
>
>
> Out of idle curiosity, I looked at Geometry Viewer.   It allows just one
> data display layer so you can only look at one table/query at a time. The
> style of the geometric objects is hard coded so you can't change their
> color, background pattern, line style etc.   It does not support color
> scales.  I can't think of any use case for which it would be useful.
>
> I would strongly suggest using QGIS instead.  Not only does QGIS have all
> of the capabilities missing from Geometry Viewer but it will read directly
> from PostGIS tables.  Since it's so widely used there's lots of
> documentation and other learning materials; I counted five different
> Youtube beginner's tutorials.
>
> Ruven Brooks
>
>
>
> On 3/8/2022 5:02 AM, Shaozhong SHI wrote:
>
> Sometime, when clicking on the eye button in PgAdmin, Geometry Viewer gets
> into a mode to display background map with options like Empty, Street,
> Topography, Grey Scale, Light Colour and Dark Matter.
>
> I am interested in making use of visualisation within PgAmin.
>
> Are there any documentation or guides on how to make the best use of these
> in PgAdmin?
>
> Regards,
>
> David
>
> ___
> postgis-users mailing 
> listpostgis-users@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/postgis-users
>
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
>
>
> ___
> postgis-users mailing 
> listpostgis-users@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/postgis-users
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Fast spatial join point to raster

[David Haynes]

Ok thanks

I aligned my query to the functional indices and it is much faster.

inner join mn_dasymetric r on ST_intersects(st_envelope(r.rast),
ST_Transform(p.geom, 26915))

On Sun, Apr 18, 2021 at 1:21 PM Paul Ramsey 
wrote:

> The indexing of rasters is (and remains a pain in the butt) (maybe
> something I should look at shortly) because the index is bound to geometry
> and almost always declraed functionally. You've built functional indexes,
>
> GIST (ST_Transform(ST_Envelope(rast),26915))
>
> But you aren't calling them with queries that use the SAME functional
> signature. The rule of functional indexes:
>
> * An index like this:  CREATE INDEX foo_x ON (function(bar))
> * must be called like this: SELECT * FROM foo WHERE function(bar) = 'baz'
>
> Note that function(bar) appears in the CREATE and in the SELECT.
>
> P
>
> > On Apr 18, 2021, at 10:53 AM, David Haynes  wrote:
> >
> > Hello,
> >
> > I'm following up on this blog post to see if there is a fast way to join
> points to rasters.
> > http://blog.cleverelephant.ca/2018/09/postgis-external-storage.html
> >
> > I'm attempting to do a dasymetric mapping process, which is assigning
> areal unit values to raster pixels. For this example, I have created two
> tables. Table 1, mn_smokers has are geom (x,y) coordinates and spatial
> indices on 4326 and 26915. The resulting raster will be stored in
> mn_dasymetric and it also has indices at 4326 and 26915. When running the
> explain query it seems that index isn't really being used.
> >
> >
> > CREATE INDEX mn_smokers_gist_4326  ON mn_smokers  USING GIST (geom);
> > CREATE INDEX mn_smokers_gist_26915  ON mn_smokers  USING GIST
> (ST_Transform(geom,26915));
> > CREATE INDEX mn_dasymetric_gist_4326 ON mn_dasymetric  USING GIST
> (ST_Transform(ST_Envelope(rast),4326));
> > CREATE INDEX mn_dasymetric_gist_26915  ON mn_dasymetric  USING GIST
> (ST_Transform(ST_Envelope(rast),26915));
> >
> > Type: Nested Loop (Inner); ; Cost: 6.22 - 1594516.11
> > Type: Limit; ; Cost: 5.11 - 6.22
> > Type: Nested Loop (Inner); ; Cost: 5.11 - 471.95
> > Type: Seq Scan; Rel: minnesota_counties ; Cost: 0.00 - 34.09
> > Type: Bitmap Heap Scan; Rel: mn_smokers ; Cost: 5.11 - 437.50
> > Type: Bitmap Index Scan; Rel: mn_smokers_gist_4326 ; Cost: 0.00 - 5.10
> > Type: CTE Scan; Rel: points ; Cost: 0.00 - 4.00
> > Type: Materialize; ; Cost: 0.00 - 59.57
> > Type: Seq Scan; Rel: mn_dasymetric ; Cost: 0.00 - 47.38
> >
> >
> > with points as
> > (
> > select sp_id as point_id, smoker as point_value, ms.geom
> > from mn_smokers ms
> > inner join minnesota_counties mc on ST_Intersects(ms.geom, mc.geom)
> > )
> > , point_raster_agg as
> > (
> > select geom, count(point_id) as pixel_point_value
> > from points
> > group by geom
> > )
> > , output_rast as
> > (
> > select r.rid, pixel_point_value,  ST_SetValue(rast, ST_Transform(p.geom,
> 26915), pixel_point_value ) as setrast
> > ,rast, geom
> > from point_raster_agg p
> > inner join mn_dasymetric r on ST_intersects(r.rast, ST_Transform(p.geom,
> 26915))
> > ___
> > postgis-users mailing list
> > postgis-users@lists.osgeo.org
> > https://lists.osgeo.org/mailman/listinfo/postgis-users
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Fast spatial join point to raster

[David Haynes]

Hello,

I'm following up on this blog post to see if there is a fast way to join
points to rasters.
http://blog.cleverelephant.ca/2018/09/postgis-external-storage.html

I'm attempting to do a dasymetric mapping process, which is assigning areal
unit values to raster pixels. For this example, I have created two tables.
Table 1, mn_smokers has are geom (x,y) coordinates and spatial indices on
4326 and 26915. The resulting raster will be stored in mn_dasymetric and it
also has indices at 4326 and 26915. When running the explain query it seems
that index isn't really being used.


CREATE INDEX mn_smokers_gist_4326  ON mn_smokers  USING GIST (geom);
CREATE INDEX mn_smokers_gist_26915  ON mn_smokers  USING GIST
(ST_Transform(geom,26915));
CREATE INDEX mn_dasymetric_gist_4326 ON mn_dasymetric  USING GIST
(ST_Transform(ST_Envelope(rast),4326));
CREATE INDEX mn_dasymetric_gist_26915  ON mn_dasymetric  USING GIST
(ST_Transform(ST_Envelope(rast),26915));

Type: Nested Loop (Inner); ; Cost: 6.22 - 1594516.11
Type: Limit; ; Cost: 5.11 - 6.22
Type: Nested Loop (Inner); ; Cost: 5.11 - 471.95
Type: Seq Scan; Rel: minnesota_counties ; Cost: 0.00 - 34.09
Type: Bitmap Heap Scan; Rel: mn_smokers ; Cost: 5.11 - 437.50
Type: Bitmap Index Scan; Rel: mn_smokers_gist_4326 ; Cost: 0.00 - 5.10
Type: CTE Scan; Rel: points ; Cost: 0.00 - 4.00
Type: Materialize; ; Cost: 0.00 - 59.57
Type: Seq Scan; Rel: mn_dasymetric ; Cost: 0.00 - 47.38


with points as
(
select sp_id as point_id, smoker as point_value, ms.geom
from mn_smokers ms
inner join minnesota_counties mc on ST_Intersects(ms.geom, mc.geom)
)
, point_raster_agg as
(
select geom, count(point_id) as pixel_point_value
from points
group by geom
)
, output_rast as
(
select r.rid, pixel_point_value,  ST_SetValue(rast, ST_Transform(p.geom,
26915), pixel_point_value ) as setrast
,rast, geom
from point_raster_agg p
inner join mn_dasymetric r on ST_intersects(r.rast, ST_Transform(p.geom,
26915))
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] issue with big data handling

[David Haynes]

Try using citusDB and sharing your data.

On Tue, Oct 29, 2019, 9:33 AM Shaozhong SHI  wrote:

> One of key issue with big data handling is that we often get the following
> message:
>
> server closed the connection unexpectedly.
>
> Any good ways to get around the problem?
>
> Regards,
>
> Shao
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] PostGIS setup on Postgres-XL

[David Haynes]

It might be easier to use citusDB.

On Wed, Jul 31, 2019, 5:09 AM Pradeep Kanth G 
wrote:

> Dear All,
>
> Recently, I have enable the Postgres-XL cluster setup but unable to
> configure the postGIS.
> Kindly help us how to configure on the cluster.
>
>
>
> Regards
>
> Pradeep Kanth
>
>
> *Disclaimer *| The information contained in this electronic message
> (including any attachments) is intended for the exclusive use of the
> addressee(s) and may contain confidential or privileged information.If you
> have received this in error, please notify the sender immediately and
> delete the material from your machine. Any  action including review,
> retransmission, dissemination of this email or the attachments present
> along with the email by persons or entities other than the intended
> recipient is prohibited.
>
>
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Error installing PostGIS on agensdb

[David Haynes]

Can you elaborate a bit more on why "the dll can not be loaded"

I'm confused as the error asks for a file at a specific location and the
file is there at that exact location.

On Wed, May 1, 2019 at 12:55 PM Regina Obe  wrote:

> Not familiar with how AgensDB puts things.
>
> That error happens whenever the dll can't be loaded.
>
>
>
> You might be better off using the binaries instead of the installer as I
> think the paths you put things relative to PostgreSQL install folder may be
> different.
>
>
>
> For example in BigSQL things go  in
>
>
>
> the lib files go into lib/postgresql folder, share/extension files in
> share/postgresql  and bin in bin
>
>
>
> So that's why the installer can never work for BigSQL and you have to use
> the binaries.
>
>
>
> The other issue is I don't ship files in the zip or installer that I know
> EDB ships to prevent accidental overwriting.
>
> So you maybe need to copy those separately from elsewhere.
>
>
>
> If you use Dependency Walker (x64 version) -
> http://www.dependencywalker.com/  will tell you the dependencies you need
> or are missing.
>
>
>
> There is yet another reason why this error happens, and that depends
> compile switches (--disable-float8-byval) comes to mind
>
> For example EDB PostgreSQL used to have this off pre-9.5 so I had (
> --disable-float8-byval  ) when compiling and so the pre-9.5 binaries I
> built weren't compatible with BigSQL  because I had to use this switch too.
>
>
>
> Post 9.5 both EDB and BigSQL don't use this switch so the binaries are
> compatible.  So if agent is doing something else, you probably can't use
> these binaries at all and will have to compile your own or pay someone to
> do so (like me J ).
>
>
>
> BTW just uploaded the 2.5.2 zip and binaries for 10.
>
>
>
> Hope that helps,
>
> Regina
>
>
>
>
>
> *From:* postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] *On
> Behalf Of *David Haynes
> *Sent:* Wednesday, May 01, 2019 12:15 PM
> *To:* postgis-users@lists.osgeo.org
> *Subject:* [postgis-users] Error installing PostGIS on agensdb
>
>
>
> Hello,
>
>
>
> I'm running into an error trying to enable the postgis extension on
> AgensDB (a fork of PostgreSQL). My setup is on a windows machine, PSQL
> 10.4. I have done this sucessfully once before, but I am unable to
> replicate it on my new machine.
>
>
>
> The process I used is the following. I downloaded the precompiled windows
> binaries (postgis-bundle-pg10x64-setup-2.5.1-1.exe) and executed it within
> the agens directory successfully, but it still errors. The error below is
> strange because the file exists within the directory.
>
>
>
> agens -p 5433 -U agens -c "CREATE EXTENSION postgis"
>
> ERROR:  could not load library
> "C:/AgensGraph-2.1.0/lib/rtpostgis-2.5.dll": The specified module could not
> be found.
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Error installing PostGIS on agensdb

[David Haynes]

Hello,

I'm running into an error trying to enable the postgis extension on AgensDB
(a fork of PostgreSQL). My setup is on a windows machine, PSQL 10.4. I have
done this sucessfully once before, but I am unable to replicate it on my
new machine.

The process I used is the following. I downloaded the precompiled windows
binaries (postgis-bundle-pg10x64-setup-2.5.1-1.exe) and executed it within
the agens directory successfully, but it still errors. The error below is
strange because the file exists within the directory.

agens -p 5433 -U agens -c "CREATE EXTENSION postgis"
ERROR:  could not load library "C:/AgensGraph-2.1.0/lib/rtpostgis-2.5.dll":
The specified module could not be found.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Improvement suggestion

[David Haynes]

Question about the solution that was posted below. I have always used
ST_CollectionExtract( ST_MakeValid(geom), 3 )  and in the post it uses
st_multi(st_makevalid(geom))

Which is preferred. On the website for ST_CollectionExtract() there is a
new warning.
When specifying 3 == POLYGON a multipolygon is returned even when the edges
are shared. This results in an invalid multipolygon for many cases such as
applying this function on an ST_Split
 result.

Try the solution outlined here:
https://gis.stackexchange.com/questions/31310/acquiring-arcgis-like-speed-in-postgis


On Sun, Dec 2, 2018 at 3:40 PM Paul Ramsey 
wrote:

> Try the solution outlined here:
>
>
> https://gis.stackexchange.com/questions/31310/acquiring-arcgis-like-speed-in-postgis
>
>
> On Sun, Dec 2, 2018 at 10:44 AM Paul van der Linden <
> paul.doskabou...@gmail.com> wrote:
>
>> As I am working with large polygons, I'm always struggling with
>> performance, and trying to find ways to improve them.
>> F.e. I have lots of queries like:
>> SELECT ST_Intersection(table1.geom,table2.geom)
>> FROM table1
>> JOIN table2 on ST_Intersects(table1.geom,table2.geom)
>>
>> In case of large polygons this is sometimes a bottleneck, and I have the
>> following suggestion:
>> Create a function which returns the relation between 2 polygons (within,
>> intersects or disjunct) so that I can do the following:
>>
>> SELECT
>>   CASE
>>  WHEN ST_Relate(table1.geom,table2.geom)=intersects THEN
>> ST_Intersection(table1.geom,table2.geom)
>>  ELSE table1.geom
>>   END
>> FROM table1
>> JOIN table2 on ST_Relate(table1.geom,table2.geom) IN (intersects,within)
>>
>> or (because ST_Relate is calculated twice in previous query):
>>
>> SELECT
>>   CASE
>>  WHEN relate=intersects THEN ST_Intersection(t1geom,t2geom)
>>  ELSE t1geom
>>   END
>> FROM (
>>   SELECT ST_Relate(table1.geom,table2.geom) as relate,table1.geom AS
>> t1geom,table2.geom AS t2geom FROM table1
>>   JOIN table2 on table1.geom && table2.geom
>> ) AS allpolies
>> WHERE relate IN (intersects,within)
>>
>> ___
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Creating postGIS extension

[David Haynes]

Hello,

I'm having a bit of trouble getting the PostGIS extension on a new VM. It
seems there were some earlier conflicts between gdal 2.0, but I have
resolved them. However, I can't create the PostGIS extension.

psql (PostgreSQL) 10.5 (Ubuntu 10.5-1.pgdg14.04+1)
david@david-VirtualBox:~$ psql -d research -c "CREATE EXTENSION postgis;"

ERROR:  could not open extension control file
"/usr/share/postgresql/9.3/extension/postgis.control": No such file or
directory

The PostGIS control file is at
/usr/share/postgresql/10/extension/postgis.control.
How can I point the server to this directory? Should I remove the older
versions of PostgreSQL and leave only version 10. There isn't any data
loaded.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] parallel spatial calculation

[David Haynes]

Thanks for the input.

My understanding of the parallel capabilities of PostgreSQL is somewhat
limited, but I believe the workers are done at the query level, without
partitioning. This likely result in different query performances compared
to a dataset that is partitioned. Also, the partitions scheme used on
spatial dataset can have a great effect on the performance of the spatial
query that you run. It is something I'm very interested in and I didn't
know if anyone in this group was working in this area.


On Thu, Aug 23, 2018 at 7:15 AM Cedric Duprez  wrote:

> Hi David,
>
> I remember two good posts on Paul Ramsey’s blog explaining query
> parallelism and PostGIS :
>
>- http://blog.cleverelephant.ca/2017/10/parallel-postgis-2.html
>- http://blog.cleverelephant.ca/2017/11/parallel-postgis-2A.html
>
> I hope it helps you.
>
> Regards.
>
> Le 22/08/2018 à 22:09, David Haynes a écrit :
>
> Hello,
>
> I have a few questions regarding parallel spatial calculations.
>
> Since PostgreSQL 10 is allowing for multiple workers is there any work on
> allowing parallel spatial calculations? Is there anyone potentially looking
> at VoltDB as a potential place for testing any parallel spatial
> calculations. VoltDB's spatial capabilities are completely undeveloped but
> it might be a place to start.
>
>
> ___
> postgis-users mailing 
> listpostgis-users@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/postgis-users
>
> ​
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] parallel spatial calculation

[David Haynes]

Hello,

I have a few questions regarding parallel spatial calculations.

Since PostgreSQL 10 is allowing for multiple workers is there any work on
allowing parallel spatial calculations? Is there anyone potentially looking
at VoltDB as a potential place for testing any parallel spatial
calculations. VoltDB's spatial capabilities are completely undeveloped but
it might be a place to start.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] raster operations

[David Haynes]

 I have a question about how to focal operations for raster statistics are
working. I started off comparing the old depreciated
function st_mapalgebrafctngb with the current version st_mapalgebra. The
difference I have determined is that the st_mapalgebrafctngb will pass the
call back function to the raster dataset similar to the st_mapalgebra, but
as the function approaches the edge of the raster tile things get wonky.
What are values that are assigned to pixels that do not have sufficient
neighbors? When I run this bit of code. The ST_Value function reports that
no values have been assigned. I assume this can be remedied by adding
information through userargs.

with smoothed_raster as
(
SELECT st_mapalgebrafctngb(rast, 1, NULL, 1, 1, 'st_mean4ma(double
precision[][][],text,text[])'::regprocedure, 'ignore', NULL) as rast
FROM glc2000_clipped_250 r
WHERE r.rid = 25
), original_raster as
(
SELECT r.rast
FROM glc2000_clipped_250 r
WHERE r.rid = 25
)
SELECT x, y, ST_Value(s.rast, x,y) as smoothed, ST_Value(o.rast, x,y) as
original
FROM smoothed_raster s, original_raster o
CROSS JOIN generate_series(1,10) as x
CROSS JOIN generate_series(1,10) as y

This leads me to a larger question regarding focal analysis functions for
raster datasets. Do these functions cross the tiles? Assume you have 9
tiles and each tile is composed of 3x3 pixels. Tiles 1-3 in row 1, 4-6 row
2 etc. If you have a neighborhood analysis function performing on tile 5.
Will it grab adjacent pixels in tile 6? Unfortunately the examples I
created returned with the same result.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Applying a formula to more than two Raster Layers

[David Haynes]

You can use the map algebra function combined with Common Table expressions.


On Tue, May 9, 2017 at 4:20 AM, David Westcott <
david.westc...@uk.rsagroup.com> wrote:

> Dear Post GIS Users
>
>
>
> Is there a way in PostGIS of applying a formula to more than two Raster
> layers?
>
>
>
> For example, I want to derive a new Raster Layer calculated as follows:
>
>
>
> Constant1 x Layer1 x (Layer2 ^ Constant2) + Constant3 x Layer3 x (Layer4 ^
> Constant4) + Constant5 x Layer5 x (Layer6 ^ Constant6)
>
>
>
> Where  Constant1, Constant2, Constant3, Constant4, Constant5 & Constant6
> are fixed constants.
>
> Layer1, Layer2, Layer3, Layer4, Layer5 & Layer6 are
> already derived Raster Layers.
>
> The ‘^’ sign represents raised to the power of. (i.e. 3^2
> = 9).
>
>
>
> How would I apply such a formula in PostGIS?
>
>
>
> Kind Regards
>
>
>
> David Westcott
>
>
>
> Please consider the environment before printing this email. RSA are proud
> to be a responsible business - www.rsagroup.com/corporate-responsibility
>
> 
> 
> *
> Royal & Sun Alliance Insurance plc (No. 93792). Registered in England &
> Wales at St. Mark's Court, Chart Way, Horsham, West Sussex, RH12 1XL.
>
> Authorised by the Prudential Regulation Authority and regulated by the
> Financial Conduct Authority and the Prudential Regulation Authority.
> For your protection, telephone calls may be recorded and monitored.
> The information in this e-mail is confidential and may be read, copied or
> used only by the intended recipients. If you have received it in error
> please contact the
> sender immediately by returning the e-mail. Please delete the e-mail and
> do not disclose any of its contents to anyone.
> No responsibility is accepted for loss or damage arising from viruses or
> changes made to this message after it was sent.
> 
> 
> *
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] PostGIS Rasters Overview Factor

[David Haynes]

You are correct overviews are powers based. Also it does not look like you
specified tile size in your raster2pgsql statement


raster2pgsql -s  -d -Y -e -I -C -F -M -t  x  -l
2,4,8,16,32

On Wed, May 10, 2017 at 11:04 AM, Osahon Oduware 
wrote:

> Hi All,
>
> I tried loading a raster with overviews using the raster2pgsql tool using
> the syntax below:
>
> raster2pgsql -s  -d -Y -e -I -C -F -M -l 2,4,8,16,32,64,128,256,512,
> 1024,2048,4096,8192,16384,32768,65536 /path/to/raster/file  |
> psql -h  -U postgres -p 5432 -d 
>
> but it gave an error message stating that the overview factor cannot be
> more than 1,000.
>
> I would like to know how the Overview-factor works. Must the value for the
> Overview-factor be in Powers of 2 (i.e. 2,4,8,16,...)?
>
> I would be glad if someone could help my understanding.
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Joining multiple rasters in a single query

[David Haynes]

Hello,

I have a question regarding the use of joins when clipping multiple raster
datasets in a single query. I have two raster datasets with the same
spatial resolution in separate tables. I perform the clips on them
separately in their own CTEs and then in the last step I would like to do
to use the map algebra to multiply the rasters together. What I am finding
is that the LEFT join is giving the correct answer and the INNER join is
not. I find this strange as both sets should be equal.

Are these assumptions correct?
The clip function should return the same number of tiles, with the same
spatial extents if the same boundary and raster datasets (spatial extent
and resolution) are used? To perform a join on raster tiles one needs to
determine if they have the same spatial extent. Currently I am using the id
fields and the ST_Within(a,b) and ST_Within(b,a) functions. (
http://postgis.net/docs/manual-1.4/ST_Equals.html)

Consider the toy queries below.
With r1 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster1 r on ST_Intersects(r.rast, p.geom)
),
r2 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster2 r on ST_Intersects(r.rast, p.geom)
),
Select r1.id, r1.name, (ST_SummaryStatsAgg(ST_MapAlgebra(r1.rast, 1,
r2.rast, 1, '[rast1]*[rast2]', '32BF'),1, True)).*
FROM r1 left join r2 on r1.id = r2.id and ST_Within(r1.rast, r2.rast)
and ST_Within(r2.rast, r1.rast)


With r1 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster1 r on ST_Intersects(r.rast, p.geom)
),
r2 as
(
SELECT p.id, p.name, ST_Clip(r.rast, p.geom) as rast
FROM polygon p inner join raster2 r on ST_Intersects(r.rast, p.geom)
),
Select r1.id, r1.name, (ST_SummaryStatsAgg(ST_MapAlgebra(r1.rast, 1,
r2.rast, 1, '[rast1]*[rast2]', '32BF'),1, True)).*
FROM r1 inner join r2 on r1.id = r2.id and ST_Within(r1.rast, r2.rast)
and ST_Within(r2.rast, r1.rast)
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] GeosIntersects Error

[David Haynes]

Hello,

I am running into an interesting version of the bad topology error.
During a raster/vector overlay I receive the following error:
"GEOSIntersects: TopologyException"

I figure I have some bad topologies in my vector dataset. I run the
following code and receive no records. I don't have any invalid topologies.

with data as
(
SELECT sgl.id as sample_geog_level_id, gi.id as geog_instance_id, gi.label
as geog_instance_label, gi.code as geog_instance_code,
ST_IsValid(bound.geog::geometry) as valid
FROM sample_geog_levels sgl
inner join geog_instances gi on sgl.id = gi.sample_geog_level_id
inner join boundaries bound on bound.geog_instance_id = gi.id
WHERE sgl.id = 1043
)
SELECT sample_geog_level_id
FROM data
WHERE valid = FALSE

However whenever I run this I receive the invalid topologies error

WITH geographic_boundaries AS
(
SELECT sgl.id as sample_geog_level_id, gi.id as geog_instance_id, gi.label
as geog_instance_label, gi.code as geog_instance_code, bound.geog::geometry
as geom
FROM sample_geog_levels sgl
inner join geog_instances gi on sgl.id = gi.sample_geog_level_id
inner join boundaries bound on bound.geog_instance_id = gi.id
WHERE sgl.id = 1043
)
SELECT p.geog_instance_id as geog_id, p.geog_instance_label as place_name,
p.geog_instance_code as place_code, ST_Clip(r.rast, 1, p.geom, True) as rast
FROM geographic_boundaries p inner join landcover.glc2000_tp  r on
ST_Intersects(r.rast, 1, p.geom)

I also went through the process of running ST_MakeValid(geom) on my data
into a new table and still encounter the error. What am I missing? This
particular error also spits out an x,y location that changes
-113.02232164290001 43.120535816300048
-126.41517874290001 47.584821516300046
-122.12053589950004 48.281250085500034
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] PL/Python

[David Haynes]

Hello Remi,

Do you know who, what list I should report this to? In my work I found is
that plpython doesn't work well with the memory drivers for gdal or ogr.
Specifically when I would run the my code it might create the raster
dataset in memory and the for no apparent reason it would terminate the
database connection. I thought this was specific to gdal memory until I
encountered some downstream issues with vector datasets.

The error was resolved choosing actual directories and writing the data
there. I suppose that there would have to a specific plpython memory driver.



On Sat, May 21, 2016 at 5:04 AM, Rémi Cura <remi.c...@gmail.com> wrote:

> Hey,
> debuggin plpython is a bit of a pain.
> Maybe you are aware that you can use
> plpy.notice to print values during execution.
>
> This blog
> <http://geoexamples.blogspot.fr/2012/01/creating-files-in-ogr-and-gdal-with.html>
> seems to indicate that what you do should work:
>
> driver = gdal.GetDriverByName( 'MEM' )
> ds = driver.Create( '', 255, 255, 1, gdal.GDT_Int32)
>
>
> Maybe you can check that executing that in plain python while logged as 
> postgres works?
>
> Cheers,
> Rémi-C
>
>
>
>
>
> 2016-05-18 21:26 GMT+02:00 David Haynes <hayne...@gmail.com>:
>
>> Yes because the second function in my doesn't error.
>>
>> On Wed, May 18, 2016 at 12:06 PM, Martijn Meijers <b.m.meij...@tudelft.nl
>> > wrote:
>>
>>> I think the gdal.getDriver call returns None. Hence the attribute error
>>> on the line afterwards. Are you sure that 19 is what you should pass as
>>> input?
>>>
>>> Martijn
>>>
>>> Verzonden vanaf mobiel.
>>>
>>>
>>>  Oorspronkelijk bericht 
>>> Van: David Haynes
>>> Datum:18-05-2016 17:35 (GMT+01:00)
>>> Aan: postgis-users@lists.osgeo.org
>>> Onderwerp: [postgis-users] PL/Python
>>>
>>> Hello,
>>>
>>> I have a question regarding importing the gdal library using pl/python.
>>> Specifically, I want to create an in memory raster using the gdal library.
>>> This is a snippet of code that I have written,
>>>
>>> CREATE OR REPLACE FUNCTION terrapop_area_level_to_raster2(
>>> sample_geog_level_id bigint, raster_variable_id bigint, raster_band bigint)
>>> RETURNS text AS
>>> $BODY$
>>> from osgeo import gdal
>>> import numpy as np
>>> from osgeo import ogr, osr
>>>
>>> memDriver = gdal.GetDriver(19)
>>> #memDriver = gdal.GetDriverByName('MEM')
>>> memRast = memDriver.Create('', 10, 10, 1, gdal.GDT_Int32)
>>>
>>> This is the error I receive. A python NoneType error. Which seems to
>>> that the gdal class has not been imported. I have verified that gdal is
>>> available on the system.
>>>
>>> ERROR: AttributeError: 'NoneType' object has no attribute 'Create'
>>> SQL state: XX000
>>>
>>> However, this function on the same server and database returns to me all
>>> the gdal drivers. Any idea how I can diagnose this problem?
>>>
>>> CREATE OR REPLACE FUNCTION gdal_drivers()
>>>   RETURNS SETOF text AS
>>> $BODY$
>>>
>>> from osgeo import gdal
>>> driver = gdal.GetDriverByName('MEM')
>>> rast = driver.Create('', 2, 4, 1, gdal.GDT_Int32)
>>>
>>> names = []
>>> driverall = gdal.GetDriverCount()
>>> for i in range(gdal.GetDriverCount()):
>>> driver = gdal.GetDriver(i)
>>> drivername = driver.ShortName
>>> names.append([i,drivername])
>>>
>>> return names
>>>
>>> $BODY$
>>>   LANGUAGE plpythonu VOLATILE
>>>
>>>
>>>
>>>
>>>
>>>
>>> ___
>>> postgis-users mailing list
>>> postgis-users@lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>>>
>>
>>
>> ___
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>>
>
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] PL/Python

[David Haynes]

Hello,

I have a question regarding importing the gdal library using pl/python.
Specifically, I want to create an in memory raster using the gdal library.
This is a snippet of code that I have written,

CREATE OR REPLACE FUNCTION terrapop_area_level_to_raster2(
sample_geog_level_id bigint, raster_variable_id bigint, raster_band bigint)
RETURNS text AS
$BODY$
from osgeo import gdal
import numpy as np
from osgeo import ogr, osr

memDriver = gdal.GetDriver(19)
#memDriver = gdal.GetDriverByName('MEM')
memRast = memDriver.Create('', 10, 10, 1, gdal.GDT_Int32)

This is the error I receive. A python NoneType error. Which seems to that
the gdal class has not been imported. I have verified that gdal is
available on the system.

ERROR: AttributeError: 'NoneType' object has no attribute 'Create'
SQL state: XX000

However, this function on the same server and database returns to me all
the gdal drivers. Any idea how I can diagnose this problem?

CREATE OR REPLACE FUNCTION gdal_drivers()
  RETURNS SETOF text AS
$BODY$

from osgeo import gdal
driver = gdal.GetDriverByName('MEM')
rast = driver.Create('', 2, 4, 1, gdal.GDT_Int32)

names = []
driverall = gdal.GetDriverCount()
for i in range(gdal.GetDriverCount()):
driver = gdal.GetDriver(i)
drivername = driver.ShortName
names.append([i,drivername])

return names

$BODY$
  LANGUAGE plpythonu VOLATILE
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Error on a raster table restore

[David Haynes]

The problem is with pg_restore for raster tables not in the public schema
use this url for reference. https://trac.osgeo.org/postgis/ticket/2485#no1

This is a fix I have implemented that seems to work, I might have all of
the raster_constraint functions though...
SHOW search_path'

Alter the following functions to include the pg_catalog in the search path.
Alter Function _raster_constraint_nodata_values (raster) SET
Search_path="$user", pg_catalog,public;
Alter Function _raster_constraint_out_db (raster) SET Search_path="$user",
pg_catalog,public;
Alter Function _raster_constraint_pixel_types(raster) SET
Search_path="$user", pg_catalog,public;
Alter Function _overview_constraint(raster, integer, name, name, name)  SET
Search_path="$user", pg_catalog,public;

now try restore command



On Wed, Jan 20, 2016 at 10:35 AM, Cedric Duprez 
wrote:

> Yes, the function ST_BandMetadata, with the right parameters, is in my
> public functions list.
> Both full versions of postgis and postgresql are the same.
>
>
> Le 15/01/2016 10:21, Tom van Tilburg a écrit :
>
> Just checking the obvious:
> Is there a function ST_BandMetedata in your functions list?
> And what is the full version output of the original server?
>
>
> On Fri, 15 Jan 2016 at 10:03 Cedric Duprez < 
> cedric.dup...@ign.fr> wrote:
>
>> Hi Tom,
>>
>> Thanks for your answer.
>> I am sure that rasters are enabled in my postgis dbase (the first thing
>> I checked). Raster functions are in the public schema.
>> This is the result of SELECT postgis_full_version() :
>> POSTGIS="2.1.0 r11822" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6
>> March 2012" GDAL="GDAL 1.11.0, released 2014/04/16" LIBXML="2.9.1"
>> LIBJSON="UNKNOWN" TOPOLOGY RASTER
>>
>> Cedric
>>
>> Le 15/01/2016 09:09, Tom van Tilburg a écrit :
>> >
>> > Cedric,
>> >
>> > Are you sure rasters are enabled in your new postgis dbase? Please
>> > check for the existence  of raster functions (st_bandmetadata in
>> > particular) in public.functions and run SELECT postgis_full_version()
>> >
>> > Tom
>> >
>> >
>> >
>> > ___
>> > postgis-users mailing list
>> > postgis-users@lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/postgis-users
>> ___
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>
>
>
> ___
> postgis-users mailing 
> listpostgis-users@lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/postgis-users
>
>
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Geoprocessing & BigData

[David Haynes]

The way that were are thinking of handling something like that now would
load the partitioned data into nodes. I imagine you want to load all the
data into 1 table on 1 node.
This would result in a large table that might be difficult for simple
queries, alternatively you could use the postgresql inheritance and load
data and create indices by state.

Address the comment that Remi-C brought up.
I have spent a lot of time investigating Postgres-Xc / Postgres-XL /
CitusDB and greenplum. None of them really parallelize your spatial
functions. This is what we have found.

The current implementation of Paragon uses a variant of round-robin
declustering. The declustering algorithm produced 1024 spatial partitions
after processing the dataset in Table 1. The physical storage and
management of the partitions in Paragon is done by taking advantage of
PostgreSQL’s sharding feature [16]. We extended the SQL create table
statement to specify spatial declustering parameters, such as, the number
of partitions to be created, the declustering method, and a label for the
declustering scheme. To execute a spatial join, the labels of the two
tables, being joined, must match. This mechanism allows the same spatial
dataset to be partitioned using different declustering schemes.

Table 1. Spatial data used for comparison

*Database Table (acronym)*

*Geometry*

*Number of Objects*

Area-water (Aw)

Polygon

39,334

Area-landmass (AI)

Polygon

5,5951

Edge (Ed)

Polyline

4,173,498



Table 2. Comparison of Query Times: Paragon vs PostgreSQL

*Query (acronym)*

*PostgreSQL (seconds)*

*Paragon (seconds)*

*Speedup*

Polygon overlaps Polygon (Aw_ov_Aw)

  77.3

  53.5

 1.37

Polyline Touches Polygon (ED_to_Al)

452.9

246.0

 1.84

Polyline Crosses Polyline (Ed_cr_Ed)

  1693.2

  1022.0

 1.65



We executed spatial join queries from the Jackpine spatial database
benchmark [6] with Paragon in a two node cluster. The queries are expressed
in SQL with some of the spatial predicates adopted by Open Geospatial
Consortium (OGC). For instance, Code 1 demonstrates the “Polyline Touches
Polygon” query shown in Table 2.

Code 1. Spatial SQL Query

SELECT COUNT(*) FROM edges ed, arealm al WHERE ST_Touches(ed.geom, al.geom);


On Wed, Jan 27, 2016 at 8:56 PM, Lars Aksel Opsahl <lars.ops...@nibio.no>
wrote:

> Hi
>
>
> We have done some testing on this using a single Postgis server.
>
>
> -layer 1 has 7924019 rows with 11 columns and about 1 billion points.
>
> -layer 2 has 1088614 rows with 20 columns.
>
> Both layers covers all of Norway.
>
>
> I do a “esri” union in a psql function and get a new table with 27852836
> rows and 30 columns with multipolygon. The size of the new table is about
> 40 GB.
>
>
> This is done in less than 3 hours (real 152m8.248s)
>
>
> I have made a psql function called func_esri.get_esri_union that I calls
> shown below.
>
>
> psql -t -q -o /tmp/vroom2.sql sl -c"drop table IF EXISTS sl_lop.r1; drop
> table IF EXISTS sl_lop.c1; select
> func_esri.get_esri_union('org_ar5arsversjon.ar5_2013_komm_flate id
> geo','org_ar.ar250_flate sl_sdeid geo', 'sl_lop.r1','sl_lop.c1',3000,false)"
>
>
> Then I take the output from this function and uses Gnu parallel to run the
> computed sqls in 20 parallel threads.
>
>
> time parallel -j 20 psql -h vroom2 -U postgres sl -c  /tmp/vroom2.sql
>
>
> This is fast Postgis server with ssd disks and a lot of memory and cpu.
>
>
> The basic idea is that I use
> https://github.com/larsop/content_balanced_grid/ to make a grid and then
> create sqls adjusted for this grid. The size of the cells varies a lot. The
> 3000 parameter in the sql function sets the limit to max 3000 bounding box
> pr cell.
>
>
> I will post the code on git hub as soon as I have time, I need to clean it
> up and make some comments first.
>
>
> We also did a small comparison with Arcgis where we ran on a small subset
> of the tables and we got result file with 186372 rows. That took about 5
> minutes with Arcgis software and 1 minute in postgres. This test was
> running on smaller postgres database server.
>
>
> Since there are different hardware for Arcgis and Postgis I will not put
> much in to this comparison, but my point is Postgis scales very good on big
> data having the right hardware and software.
>
>
> Lars
>
>
> 
> Fra: postgis-users [postgis-users-boun...@lists.osgeo.org] på vegne av
> Ravi Pavuluri [ravith...@ymail.com]
> Sendt: 27. januar 2016 21:31
> Til: PostGIS Users Discussion
> Emne: Re: [postgis-users] Geoprocessing & BigData
>
> Hi David,
>
> I are dealing with census blocks/census block groups spanning a few
> million records.
>
> Thanks,
> Ravi.
>
> On Monday, January 25, 2016 

Re: [postgis-users] Geoprocessing & BigData

[David Haynes]

We have done some work, implementing parallel spatial queries using a
spatial declustering algorithm. How large are your datasets?

On Mon, Jan 18, 2016 at 1:51 PM, Rémi Cura  wrote:

> Hey,
> if you have one beefy server you can parallelize throwing several queries
> working on sub set of your data.
> (aka parallel processing trough data partition).
> One conceptual example : you want to process the world, you create 20
> workers, a list of countries, and then make the worker process the list
> country by country.
>
> If you think one postgres server will not be sufficient,
> you could of course shard your data across several servers,
> with options ranging from writting from scratch (you rewrite everything),
> to using existing open source code, to dedicated solution like
>  Postgresql-Xc, greenplum, ...
>
> However, sorry to say this but in your case it looks like your first
> improvement step will not come from massive paralleling but from first
> better understanding the world of geospatial data and postgis.
>
> Cheers,
> Rémi-C
>
> 2016-01-18 19:30 GMT+01:00 Vincent Picavet (ml) :
>
>> Hi Ravi,
>>
>>
>>
>>
>> On 18/01/2016 19:14, Ravi Pavuluri wrote:
>> > Hi All,
>> >
>> > I am checking if there is a way to process quickly large datasets such
>> > as census blocks in PostGIS and also by leveraging big data platform. I
>> > have few questions in this regard.
>> >
>> > 1) When I try intersect for sample census blocks with another polygon
>> > layer, PostGIS 2.2(on Postgres 9.4) takes ~60 minutes (after optimizing
>> > from http://postgis.net/2014/03/14/tip_intersection_faster/ ) while on
>> > ESRI ArcMap takes ~10 minutes. PostGIS layers already have geospatial
>> > indices. Is there anyway to optimize this further?
>>
>> Following the links on your page, here is a good answer from Paul (TL;DR
>> : st_intersection is slow, avoid it) :
>>
>> http://gis.stackexchange.com/questions/31310/acquiring-arcgis-like-speed-in-postgis/31562
>>
>> > 2) What is an equivalent of ESRI Union in PostGIS? I didn't see any out
>> > of the box functions and any tips here are appreciated.
>>
>> If ESRI Union makes a union, maybe st_union ? But I guess there are some
>> semantic issues here.
>>
>> > 3) Is there anyway we can expedite these geoprocessing
>> > tasks(union/intersect etc) using big data platform (Ex: hadoop)? Most
>> > examples talk about analysis (contains etc)  but not about geoprocessing
>> > on geospatial data. Any input is appreciated.
>>
>> Lots of people do geoprocessing too with PostGIS, including long-running
>> jobs on large volumes of data ( worldwide osm data processing namely).
>> "Big data" is a really subjective word. Are your geoprocessing needs
>> really parallelizable ? What kind of volumes are we talking about ? MB,
>> GB, TB ? What kind of hardware do you have at hand ?
>>
>> One way to do some sort of map-reduce with PostGIS is to use a bunch of
>> servers with FDW connections between a source master and these slaves,
>> map the data processing to the slave servers and reduce it on the main
>> server. With a bit of Python as glue code this can be automated and
>> quite efficient, even though this kind of sharding is not automated (
>> yet ?).
>>
>> Vincent
>>
>> >
>> > Thanks,
>> > Ravi.
>> >
>> >
>> > ___
>> > postgis-users mailing list
>> > postgis-users@lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/postgis-users
>> >
>>
>> ___
>> postgis-users mailing list
>> postgis-users@lists.osgeo.org
>> http://lists.osgeo.org/mailman/listinfo/postgis-users
>
>
>
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users
>
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Raster calculation on 3d raster (stack)

[David Haynes]

Hello,

I have a mulibanded raster and I am wondering if I can apply any of the
raster processing functions (min, max, mean, stdev) to multibanded rasters.
Specifically, I am wondering if these functions will be able to applied to
the z-axis.

Here is an example, I have worked up in python starting with two 2d arrays.
x = [[1,4,5,6], [3,5,6,2]]
z = [[2,7,8,9], [8,3,5,9]]
convert to np.array
y = np.dstack((np.array(x), np.array(z)))
>>> y.sum(2)
array([[ 3, 11, 13, 15],
   [11,  8, 11, 11]])

I want to apply the raster processing functions to the z-axis and I am
wondering if I can do it in PostGIS? It seems like the mapalgebra call back
functions could do this, but I am having difficulty finding an example.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Parallelisation provides powerful postgis performance perks (script + ppt slides) [x-posted: pgsql-performance]

[David Haynes]

Hello,

I am a researcher at the University of Minnesota, US. who is working on a
project that uses PostGIS as our platform. I have been researching a number
of PostgreSQL platforms that claim to support parallelizing PostGIS or
geographic analysis. I am interested to learn more about the wrapper you
have written. Have you looked into the other platforms such at pg_shard
(CitusDB), or postgres-xl? We have not found them effective for actually
distributing a spatial query. However, I will test your code out in our use
case.

When I briefly look at the text you have written in the Quick Example It
seems that you are distributing your query by an ID field. I am wondering
how your method would apply to raster datasets? Distributing geographic
data by an ID can get you into problems because of the dependency for
certain analytical functions.

This sounds great, hope to hear back from you soon.

On Thu, Jul 23, 2015 at 8:54 AM, Graeme B. Bell graeme.b...@nibio.no
wrote:

 Hi all,

 Do you run map intersections between national scale geometry maps in
 postgis? (or other long-running GIS operations?)
 Do you hate that feeling of waiting all day/week for the query to complete?
 Well, here's a solution for you.

 1. For those that don't like par_psql (http://github.com/gbb/par_psql),
 this alternative approach uses the Gnu Parallel command to organise
 parallelism for queries that take days to run usually. We saw up to 20x
 performance improvements here, a day's work in one hour. May give you a few
 ideas about how to parallelise your own code with Gnu Parallel.

 https://github.com/gbb/fast_map_intersection


 2. Also, I gave a talk at FOSS4G Como about these tools (and a few
 others), and how to get better performance from your PostGIS database with
 parallelisation.
 This may be helpful to people who are new to parallelisation / multi-core
 work with postgres for GIS work.

 http://graemebell.net/foss4gcomo.pdf

 Enjoy,

 Graeme Bell.
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[postgis-users] ST_Clip

[David Haynes]

Hello

I believe I am in appropriately using the ST_Clip function on a multiband
raster. I have a raster with 13 bands and I want to retrieve from the clip
bands # 1,3

How do I do that?

I tried using:
array[1,3] as bnd

Didn't work still returns 13 bands.

with bnd_num as
(
select array[1,3] as bnd
)
select p.geoid, p.label, ST_Numbands(ST_CLIP(r.rast,p.geom,b.bnd, True))
from bnd_num b, polygon p inner join modis_igbp_stack r on
ST_Intersects(r.rast, p.geom);
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] polygons/geometry files for European provinces?

[David Haynes]

A reason for a shameless plug  terrapop.org
All of its free

On Fri, Apr 10, 2015 at 1:46 PM, Martijn Meijers b.m.meij...@tudelft.nl
wrote:

 Not free:

 http://www.eurogeographics.org/products-and-services/euroboundarymap


 Openstreetmap may have what you need as well:

 http://www.itoworld.com/map/2?lon=9.60125lat=51.11984zoom=7


 Martijn

 On 10-04-15 20:07, Joseph Spenner wrote:

 Does anyone have a good source of polygon shapes which represent European
 countries/provinces?

 I found a few sites, such as:

 https://www.sciencebase.gov/catalog/items?q=filter=tags%
 3DGeologic+province

 But it does not provide the political boundaries.

 I'm looking for the provinces which can be displayed at a site such as:
 http://www.meteoalarm.eu/
   (where clicking on the counties displays the provinces; those provinces
 are what is of interest).

 Any help/links would be great.

 Thanks!

 Regards,
 Joseph Spenner






 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users


 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Fwd: Partitionning using geometry

[David Haynes]

Actually we are in the process of doing something similar

On Wed, Apr 1, 2015 at 11:03 AM, Rémi Cura remi.c...@gmail.com wrote:

 (cross-post from postgres list)
 Hey dear list,

 I'd like to partition geographical (geometry) data with postgres mechanism.
 (my usage is in fact related to pointcloud, but I use geometry as a work
 around)
 From example I read on constraint, nothing should prevent it from working
 Here is a self contained example, the planner doesnt seems to use the
 constraint_exclusion mechanism, whatever the constraint

 Thanks,
 Cheers,
 Rémi-C

 --

 CREATE SCHEMA IF NOT EXISTS test_partitionning;
 SET search_path TO test_partitionning, public ;

 DROP TABLE IF  EXISTS test_father CASCADE;
 CREATE TABLE test_father  (
 gid SERIAL PRIMARY KEY
 , geom geometry
 );

 create table test_child_1 (
 check (geometry_overlaps(geom,ST_Expand(ST_MakePoint(10,10),10  ) ) )
 ,check ( geomST_Expand(ST_MakePoint(10,10),10  ) )
 , CHECK (ST_X(geom) BETWEEN 0 AND 20)
 , CHECK (ST_Y(geom) BETWEEN 0 AND 20)
 , CHECK (  ST_Intersects(geom, ST_Expand(ST_MakePoint(10,10),10  ))  )
 ) inherits (test_father);
 --CREATE INDEX ON test_child_1 USING GIST(geom);

 create table test_child_2 (
 check (geometry_overlaps(geom,ST_Expand(ST_MakePoint(30,10),10  ) ) )
 ,check ( geomST_Expand(ST_MakePoint(30,10),10  ) )
 , CHECK (ST_X(geom) BETWEEN 20 AND 40)
 , CHECK (ST_Y(geom) BETWEEN 0 AND 20)
 , CHECK (  ST_Intersects(geom, ST_Expand(ST_MakePoint(30,10),10  ))  )
 ) inherits (test_father);
 --CREATE INDEX ON test_child_2 USING GIST(geom);


 INSERT INTO test_child_1 (geom)
 SELECT ST_MakePoint(s1/10.0+random(),s2/10.0+random())
 FROM generate_series(1,90) AS s1, generate_series(1,90) AS s2;

 INSERT INTO test_child_2 (geom)
 SELECT ST_MakePoint(s1/10.0+random(),s2/10.0+random())
 FROM generate_series(200,300) AS s1, generate_series(1,90) AS s2;


 SHOW constraint_exclusion;
 SET constraint_exclusion TO partition;


 WITH area_of_interest AS (
 SELECT ST_Buffer(ST_MakePoint(5,5),1) as buf
 )
 SELECT *
 FROM area_of_interest, test_father
 WHERE  -- geom  buf
 ST_X(geom) BETWEEN ST_XMin(buf) AND ST_Xmax(buf)
 AND ST_Y(geom) BETWEEN ST_YMin(buf) AND ST_Ymax(buf) ;


 SELECT *
 FROM  test_father , ST_Buffer(ST_MakePoint(5,5),1) as buf
 WHERE
 ST_X(geom) BETWEEN ST_XMin(buf) AND ST_Xmax(buf)
 AND ST_Y(geom) BETWEEN ST_YMin(buf) AND ST_Ymax(buf);
 --


 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] calculating area

[David Haynes]

There are many ways to determine area in PostGIS, but I will simplify this
to 3 groups.
1 consider the earth as an ellipsoid, where the radius changes, 2 consider
the world to a perfect sphere, 3 your area is small enough to fit in a
defined equal area projection. For work, I have generated a number of area
datasets for country boundries using the functions in PostGIS. You will
notice that the 3rd function uses a global equal area projection. They are
all different and in someways they are all accurate. The first function
uses a mathematical model that tries to model the earth and takes a bit
longer. The second uses a sphere and is the fastest, while the third
flattens the earth and its accuracy is based on the projection.

1) Spherical Area ST_Area(p.geog, false) as ellipsoid / spheroid
The current World Geodetic System
http://en.wikipedia.org/wiki/World_Geodetic_System model uses a spheroid
whose radius is 6,378.137 km at the equator
http://en.wikipedia.org/wiki/Equator and 6,356.752 km at the poles
http://en.wikipedia.org/wiki/Geographical_pole.
2) Spheroid AREA ST_Area(p.geog) as sphere,
3) Projected Area ST_Area(ST_Transform(p.geog::geometry,3410)) proj_area

In general, I would recommend the function 1 as an accurate calculation. I
have document of country areas I can share with you if you like. It shows a
collection of country polygons that have the area determined using all
three different methods. In this document, I am assuming that spheroid will
be the truth.
The interesting thing is that there are major differences based on the
latitude. The area of a geography goes from over to under projection at the
35 degree. This is due to shrinking and expansion of the ellipsoid and
equal area projection.

On Wed, Feb 25, 2015 at 6:53 AM, Lelo - Luiz Rogério De Pieri 
lelo.pi...@gmail.com wrote:

 Hi,

 I think you always have to transforme to UTM to have area in square meters.

 Best Regards,

 2015-02-25 9:10 GMT-03:00 Ahmet Temiz ahmettemi...@gmail.com:

 Hello,

 How  do we calculate the area (polygon) in country ( Turkey) level ?

 Which one does give better result ? transforming to UTM or using
 geography ?

 regards

 --
 Ahmet Temiz
 Jeoloji Müh.
 Afet ve Acil Durum Yönetimi Başkanlığı
 Planlama ve Zarar Azaltma Dairesi Başkanlığı


 

 Ahmet Temiz
 Geological Eng.
 Information Systems - GIS Group
 Disaster and Emergency Management
 of Presidency

 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users




 --
 Rogério De Pieri  (Lelo)
 MBA em Gerenciamento de Projetos - FGV
 SCJP 5
 Buscando melhorar a cada dia
 Áudio, Hardware  Software
 www.twitter.com/lelopieri
 blogdolelo.wordpress.com

 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] [Raster] Finding pixels intersecting an extent (polygon)

[David Haynes]

select ST_Clip(r.rast,p.geom) as rast
from polygon p inner join raster r on ST_intersects(r.rast, p.geom)

This returns a raster which has all pixels inside the polygon

On Thu, Feb 5, 2015 at 4:05 PM, Jean Marchal jean.d.marc...@gmail.com
wrote:

 Hi list,

 I am trying to return all the pixels in a raster that intersect (not just
 touch) an extent (say a rectangle). I tried ST_Clip and
 ST_Intersection(raster, geom) but they don't return all the pixels that
 intersect my extent polygon. Do I have to vectorize the raster first using
 ST_PixelAsPolygons or there is a better / more efficient way to proceed?

 Ultimately the goal is to fetch the resulting raster in R.

 Thanks,

 Jean

 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Replicating ArcGIS relationship classes in PostGIS

[David Haynes]

The only item that immediately springs to mind would be creating views to
preserve your join relationships. Any user that queried the view would
automatically have that underlying joined information available.

Perhaps you could explain more about the relationship classes you are
trying to implement.

On Thu, Oct 23, 2014 at 11:32 PM, Jochen Albrecht jochen.albre...@gmail.com
 wrote:


 I am working on behalf of a small non-profit that has so far used ArcGIS
 but is interested in moving open source. What kept them from making the
 move so far is that they are heavy users of relationship classes. The only
 discussion (inconclusive) that I could find was initiated by Lee
 Hachadoorian in 2008. Lots has changed since and I am wondering whether
 anybody has actually done this (implementing relationship classes in
 PostGIS). For example, back then, Lee suggested a trial-and-error approach
 to understand the encoding of relationship class tables. Have we gotten any
 further on this front? I am surprised that GIS Stack Exchange has no
 mentioning of this at all (assuming the search tools work).
 We are using relationship classes for 1:many and many:many relationships.
 This is not a problem for PostGIS; I just don't know how to migrate the
 geodatabase keeping those relationships intact.
 Cheers,
  Jochen

 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] PostGIS supported by which ArcGIS

[David Haynes]

I have not been able to use Arc 10.3 for visualizing raster data stored in
postgresql

On Thu, Oct 2, 2014 at 2:30 AM, Birgit Laggner birgit.lagg...@ti.bund.de
wrote:

 We use ArcGIS 10.2 for read/write interaction with our PostgreSQL
 database. We only have ArcGIS Desktop without SDE. There is a driver
 (database client file) from ArcGIS for PostgreSQL, you have to install. We
 are more or less content with this concept. Surely, the performance is not
 as if you would be working directly on the PostgreSQL database, but still,
 it isn't necessary to always export and import into and from shapefile.
 That's a huge plus in my opinion.

 Here is a link to a documentation site of ArcGIS regarding the
 requirements for the interoperability with PostgreSQL:

 http://resources.arcgis.com/en/help/system-requirements/10.2/index.html#//
 0151007500

 Regards,

 Birgit.


 Am 27.09.2014 20:28, schrieb Stefan Keller:

  Actually yes, I meant the PostGIS option of storing data in
 PostgreSQL/PostGIS from ArcGIS.

 --S.

 2014-09-27 20:20 GMT+02:00 Randal Hale rjh...@northrivergeographic.com:

 as far as I understand it (and I could be wrong on some of this)

 Desktop will support reading if you load the postgresql library files
 (that
 can be downloaded form your ESRI customer care portal).
 Once those are loaded you can read through a definition query data from
 postgis into arcgis. I just made a map using arcgis and reading data from
 postgis.
 It doesn't support write unless you have SDE for Postgresql - which isn't
 postgis. That dives into the world of ArcGIS Server


 ArcGIS supposedly is moving to being able to read/write the geopackage
 format which to me opens read/write to PostGIS. BUT if they open read
 write
 for postgis that starts a slow death of SDEso I don't believe it ever
 will (for the foreseeable future).

 Randy




 On 09/27/2014 02:10 PM, Stefan Keller wrote:

 I'm asked from time to time for an advice which ArcGIS product
 supports PostGIS (read/write)?
 Even after reading the Esri pages I'm not sure if it's only ArcGIS
 for Server Workgroup (formerly ArcSDE).
 Does it work also with ArcGIS for Desktop Basic (ArcView)?

 Yours, S.
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users


 --
 -
 Randal Hale
 North River Geographic Systems, Inc
 http://www.northrivergeographic.com
 423.653.3611 rjh...@northrivergeographic.com
 twitter:rjhale http://about.me/rjhale
 http://www.northrivergeographic.com/spatial-connect

 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users


 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[postgis-users] ST_AddBand

[David Haynes]

I have a question about this function, regarding the fifth variant.
Why does the function assume that all of the rasters input in the array
have their information stored at the same band? The variant will only allow
for 1 band number to be specified for all rasters.

raster *ST_AddBand*(raster torast, raster[] fromrasts, integer fromband=1,
integer torastindex=at_end);

This is will not work

ST_AddBand(ST_MakeEmptyraster(layers.stack_rast),ARRAY[layers.stack_rast,
layers.stack_rast], ARRAY[1,3]) as rast

This does, but assumes that I want band 1 both times.
ST_AddBand(ST_MakeEmptyraster(layers.stack_rast),ARRAY[layers.stack_rast,
layers.stack_rast], 1) as rast

Shouldn't there be a variant that allows that uses an array integer that
allows for you to individual specifications of bands.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Accessing NWS poly data with lat/lon input

[David Haynes]

Difficult to understand what output you want.

ST_AsText(geom) - will give you the bounding information about your county
polygons.
ST_AsCentroid(geom) - will give you the centroid of the polygon.

What are the lat/long values of?


On Fri, Aug 8, 2014 at 11:24 AM, Nicolas Ribot nicolas.ri...@gmail.com
wrote:

 Hello

 Considering your geometric column is named geom:

 select geom
 from polys
 where st_contains(geom, st_setSRID(st_makePoint(lon, lat), 4326));

 (if several polygons contain this point, they will all be returned. Adding
 LIMIT 1 at the end of the query will force to get only one result.

 Nicolas


 On 8 August 2014 17:24, Newcomb, Doug doug_newc...@fws.gov wrote:

 Joseph,
 I don't see a geometry or geography field in your table listing.

 Doug


 On Fri, Aug 8, 2014 at 11:12 AM, Joseph Spenner joseph85...@yahoo.com
 wrote:

 Hello, I've just got postgresql93-server installed, with
 postgis21-postgresql93 on a CentOS system.

 I then used shp2pgsql to load in some shape files using this procedure:


 http://suite.opengeo.org/docs/latest/dataadmin/pgGettingStarted/shp2pgsql.html

 The shapefiles were obtained from NWS:
 http://www.nws.noaa.gov/geodata/catalog/wsom/html/pubzone.htm

 Everything appears to have been written to the database properly:

 polydata=# select state,time_zone,zone,name,lon,lat from polys where
 state='CO' limit 10;
  state | time_zone | zone |   name   |
 lon   |  lat

 ---+---+--+--++---
  CO| M | 048  | Logan County |
 -103.110114271 | 40.7246902558
  CO| M | 044  | Morgan County|
 -103.809823690 | 40.2627093692
  CO| M | 050  | Sedgwick County  |
 -102.351810279 | 40.8758426210
  CO| M | 049  | Washington County|
 -103.201287262 | 39.9710250432
  CO| M | 090  | Yuma County  |
 -102.424258955 | 40.0029195574
  CO| M | 099  | Springfield Vicinity/Baca County |
 -102.560453567 | 37.3192132374
  CO| M | 097  | Las Animas Vicinity/Bent County  |
 -103.071690129 | 37.9551177349
  CO| M | 092  | Cheyenne County  |
 -102.603398004 | 38.8279341557
  CO| M | 089  | Crowley County   |
 -103.784878922 | 38.3266440954
  CO| M | 091  | Kit Carson County|
 -102.602884309 | 39.3054124576
 (10 rows)

 polydata=#

 My goal is to construct a query providing lat/lon such that the result
 returned is the single row poly containing that point.

 Does anyone know how I would construct this query?

 Thanks!

 Regards,
 Joseph Spenner


 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users




 --
 Doug Newcomb
 USFWS
 Raleigh, NC
 919-856-4520 ext. 14 doug_newc...@fws.gov

 -
 The opinions I express are my own and are not representative of the
 official policy of the U.S.Fish and Wildlife Service or Dept. of the
 Interior.   Life is too short for undocumented, proprietary data formats.

 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users



 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Sorry for simple question!

[David Haynes]

You can place rasters into a single raster using ST_BAND. If you are
familiar with rasters this will make a multispectral raster. I don't
believe there is a limit to the bands as I have added over 50. To make the
multispectral raster your other rasters (matrices) will need to be the same
cell size and alignment.

You might try http://www.gdal.org/gdal_merge.html (using the -separate)

This will create a 2 banded raster in postgis
ST_AddBand(ST_MakeEmptyraster(layers.modis_rast),ARRAY[layers.modis_rast,
layers.stack_rast]) as rast


On Fri, Aug 1, 2014 at 6:42 AM, Darrel Maddy darrel.ma...@newcastle.ac.uk
wrote:

  Dear all,



 I am just reacquainting myself with postgres and postgis and trying to
 decide how best to store the output from my model (large matrices stored in
 gtiff format).  I really want to store multiple rasters in one table but
 when I try to visualise these in QGIS (via an sql extract in dbmanager
 using import to new layer) nothing happens. I can add one raster per table
 to the canvas but I do not want thousands of tables.   Will the dbmanager
 sql add layer work with rast as a geometry?



 I suspect this is simple, but it eludes me at present.



 Thanks in anticipation of your help.



 Darrel

 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Reference vrt datasets in postgis raster

[David Haynes]

Hello,

You might want to try this projection:
http://spatialreference.org/ref/sr-org/6842/postgis/
This the most accurate projections, despite the description on the spatial
reference site. I am doing related work and found this calculator at the
following url
http://landweb.nascom.nasa.gov/cgi-bin/developer/tilemap.cgi
There is no projection accuracy error when using SRID 6842.



On Fri, Jul 11, 2014 at 9:33 AM, Christian Yrrman yrr...@gmail.com wrote:

 Hi All,


 I'm a newbie to this group, so hello everyone and thanks in advance for
 letting me participate :-)

 After several tests and searching on the net I'm a little stuck with the
 following problem:

 I have downloaded a bunch of MODIS tiles (MOD13A1) from [0] that I want to
 (only) reference in a postgis raster database. These tiles come in hdf
 containers with several layers in it
 .
 - For a start,  I'm only interested in the EVI layer, so I created EVI vrt
 files with command [1]
 - This gives me several EVI*.vrt files referencing EVI with the original
 sinusoidal projection. According to spatialreference.org, the closest SRS
 match is sr-org6974 [2].
 - As this is not in postgis, I've imported this into spatial_ref_sys [3]
 - Next, I'm using the following command to insert the data into raster
 table modistest with the following command:

 for i in `find / -name EVI_origSRS*.vrt`; do
 raster2pgsql -s 96974 -R -F -f rast -a $i modistest;
 done

 where:
 -s 96974 - use the freshly imported MODIS SRS from [2]
   -R  Register the raster as an out-of-db (filesystem) raster.  Provided
   raster should have absolute path to the file
   -F  Add a column with the filename of the raster.
   -f column Specify the name of the raster column
  -a  Appends raster into current table, must be
  exactly the same table schema.

 I've built the index on rast in a separate step (not using switch -I)

 This populates table modistest (using | psql...), so it looks fine.
 However when I want to retrieve data from the raster, e.g. via

 SELECT ST_AsTiff(rast) from public.modistest where id=1158;

 I get an error message like this
 ERROR:  rt_raster_from_gdal_dataset: Unable to get data from transformed
 raster
 CONTEXT:  PL/pgSQL function st_astiff line 20 at RETURN

 and in the postgresql log I get

 ERROR 4:
 `HDF4_EOS:EOS_GRID:MOD13A1.A2012065.h35v10.005.2012082112322.hdf:MODIS_Grid_16DAY_500m_VI:500m
 16 days EVI' does not exist in the file system,
 and is not recognised as a supported dataset name.


 So it looks like the layers were not correctly referenced - however I
 cannot see a mistake. I'm actually using find / -name EVI_origSRS*.vrt
 to get the full path for the loop (not for i in *vrt)


 Any help/advice is highly appreciated! In a next step, I would even like
 to transform the sinusoidal tile afterwards with ST_transform to epsg:4326.
 That would be fantastic, I'd have all tiles in original format in the file
 system and let postgis do all the rest.. :-)



 Thanks in advance and greetings!



 Chris





 [0] MODIS download website
 http://e4ftl01.cr.usgs.gov/MOLT/MOD13A1.005/

 [1] command to create vrt files:
 #!/bin/bash

 for i in *hdf; do
 echo Processing $i
 j=`gdalinfo $i| grep SUBDATASET_2_NAME | cut -d '=' -f 2`
 echo Extracted $j
 gdal_translate -of VRT $j EVI_origSRS_$i.vrt

 done

 [2] SRS projeciton for MODIS
 http://spatialreference.org/ref/sr-org/6974/

 [3] SQL command to insert MODIS SRS

 INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext)
 values ( 96974, 'sr-org', 6974, '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 
 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ',

 'PROJCS[MODIS Sinusoidal,
 GEOGCS[WGS 84,DATUM[WGS_1984,SPHEROID[WGS 
 84,6378137,298.257223563,AUTHORITY[EPSG,7030]],
 AUTHORITY[EPSG,6326]],PRIMEM[Greenwich,0,AUTHORITY[EPSG,8901]],

 UNIT[degree,0.01745329251994328,AUTHORITY[EPSG,9122]],AUTHORITY[EPSG,4326]],
 PROJECTION[Sinusoidal],
 PARAMETER[false_easting,0.0],

 PARAMETER[false_northing,0.0],
 PARAMETER[central_meridian,0.0],
 PARAMETER[semi_major,6371007.181],
 PARAMETER[semi_minor,6371007.181],
 UNIT[m,1.0],AUTHORITY[SR-ORG,6974]]');




 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[postgis-users] ST_SelectByValue

[David Haynes]

Hello,

I am fairly new to using rasters in postgis, but I am wondering if there is
a faster way to extract values from categorical rasters (e.g. landcover
types) than using map algebra. In the PostGIS in Action book, pg 404 has a
function called ST_SelectByValue(raster|geometry, 'expression') which seems
to incorporate map algebra, however I can't seem to tell if this function
is finished and available in postgis as I get an unknown function error.

My version is this.
PostgreSQL 9.3.2 on x86_64-unknown-linux-gnu, compiled by gcc
(Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit POSTGIS=2.1.1 r12113
GEOS=3.3.3-CAPI-1.7.4 PROJ=Rel. 4.7.1, 23 September 2009 GDAL=GDAL
1.9.0, released 2011/12/29 LIBXML=2.7.8 LIBJSON=UN (...)

And I found a webpage
http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking03
search for ST_SelectByValue at the bottom, indicating that it was
finished.

AAccomplished Objectives
http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking03#AccomplishedObjectives

*ST_ValueCount(raster, value) - integer - done see below*

*ST_ValuePercent(raster, value) - double precision - done see below*

*ST_Resample(raster, method, originx, originy, pixelsizex, pixelsizey) -
raster - done see below*

*ST_SelectByValue(raster, 'expression') - same type as first argument*

*ST_Clip(raster|geometry,geometry) - same type as first argument*
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users