Re: [postgis-users] Dted elevations
We solved the issue. It turns out that searching for elevation on DTED level 2
takes a long time in comparison to DTED Level 1. This is not surprising
considering the density of the data. I'm wondering if anyone has looked into
optimizing this kind of elevation searching?
--Jack
-Original Message-
From: postgis-users-boun...@lists.osgeo.org
[mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of Gold, Jack L (US
SSA)
Sent: Friday, May 09, 2014 11:01 AM
To: 'postgis-users@lists.osgeo.org'
Subject: [postgis-users] Dted elevations
I know this question has already been answered a hundred times but I am at a
remote location currently with only email access right now and I hope someone
can help. The following query is taking close to a second to run and I think
it should be much faster.
Select st_value from (with f as ( select st_transform ( st_SetSRID(
st_MakePoint('106','32'), 4326) as geom) select st_value(dted_elevations.rast,
f.geom) from dted_elevations cross join f where st_intersects (
dted_elevations.rast , f.from) and st_value (dted_elevations.rast, f.geom) is
not null) as st_value;
Any help is greatly appreciated.
Jack Gold
___
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] ST_MakeLine() woes
Sorry it didn't work :-( . Can you check that temp_id is not always 1?
It is more that this approach is flawed
_depending on the raster resolution and the small turn in your line (fixe =
simplify line).
_anyway, you are constructing a new line from centroid, and not putting the
height information on the old line (see end of the message)).
Can you try to spatially filter your line : with the same code, change just
the beginning
WITH line AS
-- From an arbitrary line
(SELECT ST_Simplify(geom,2*your_raster_cell_size) from foot_cl where
(sheet='AT24' AND sid=463)),
cells AS
A safer way to do it would be for instance:
find the pixel under the line
transform those pixels into polygons along with height value
(ST_PixelAsPolygons)
cut the line with those polygons (create new points at each enter/exit of
pixels)
add the height information of each polygon into the line parts.
renode the line by assembling in the correct order the line parts.
A suggestion of code for it:
-
WITH rast AS ( --getting a raster, and a diag line for this raster
SELECT * --,
ST_Intersection(rast,line)--unnest(regexp_split_to_array(ST_Summary(rast),'\n'))
FROM
test_raster.test_temp_raster,ST_SetSRID(ST_GeomFromText('Linestring(651050.
6860677,651055 6860682)') ,931008) AS line
WHERE rid = 273143
LIMIT 1
)
,pix_under_line AS ( --get the pixels that are covered by the line,
transform the pixels into square (pylgon), keep the value of the band 1
(heigth is supposed to be here)
--this is suboptimal and could be replacer by using
st_intersection
SELECT pix.*, line
FROM rast,ST_PixelAsPolygons(rast,1) AS pix
WHERE ST_Intersects(pix.geom,line)=TRUE
--AND pix.val!=0 --if raster from interpolation, no need to
keep wrong parts
)
,cutting_line_with_pix AS (--splitting the line with the pixels in
order to obtain multiple parts of the line, each covering one pixel
--(see
https://github.com/Remi-C/_utilities/tree/master/postgis for
rc_split_multi)
SELECT ST_Dump(rc_Split_multi(min(line),
ST_Union(ST_Boundary(geom)),0.01)) AS splitted_line
FROM pix_under_line
)
,splitted_line AS ( --filtering the obtained parts of line to remove
ghost created by precision error and point-line
SELECT (splitted_line).path, (splitted_line).geom, ST_AsText(
(splitted_line).geom)
FROM cutting_line_with_pix
WHERE ST_Length((splitted_line).geom)0.001
)
,sl_and_pix AS ( --for each parts of line, get the pixel polygon that
it covers, along with the value of the pixel. Add this value to the parts
of line
SELECT DISTINCT ON (path) sl.*, pul.*,
ST_AddMeasure(sl.geom,pul.val,pul.val) AS l_heigth
FROM splitted_line AS sl, pix_under_line AS pul
WHERE ST_Intersects(sl.geom,pul.geom)=TRUE
ORDER BY path ASC, ST_Length(sl.geom) DESC
)--fusion the line parts to create a single line (don't use ST_Union,
it drops the M value)
SELECT ST_Astext(ST_MakeLine(l_heigth ORDER BY path ASC) )
FROM sl_and_pix
-
Cheers,
Rémi-C
2014-05-13 1:08 GMT+02:00 georgew gws...@hotmail.com:
thanks Remi, attached is the result after running your code, no change
unfortunately.
I had to make some slight changes to the code to make it work. The image
has
the code.
However you also say that my approach is not as robust as it could be. Any
suggestions on how to make it more robust?
newline.png http://postgis.17.x6.nabble.com/file/n5006277/newline.png
--
View this message in context:
http://postgis.17.x6.nabble.com/ST-MakeLine-woes-tp5006266p5006277.html
Sent from the PostGIS - User mailing list archive at Nabble.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] Postgis raster tutorial
Done! http://geospatialelucubrations.blogspot.ca/2014/05/a-guide-to-rasterization-of-vector.html Pierre -Original Message- From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users- boun...@lists.osgeo.org] On Behalf Of georgew Sent: Monday, April 21, 2014 6:37 PM To: postgis-users@lists.osgeo.org Subject: Re: [postgis-users] Postgis raster tutorial Thank you Pierre, your add-ons are a life saver. I'll make sure to make the most of them. And I look forward to your blog in 2000 minutes, the countdown has started already! -- View this message in context: http://postgis.17.x6.nabble.com/Postgis- raster-tutorial-tp5006124p5006169.html Sent from the PostGIS - User mailing list archive at Nabble.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] getting raster values for point from st_value
Create a new table with the extents of the raster tiles and check if it overlaps with your points: CREATE TABLE rastextets AS SELECT rast::geometry FROM tbl_raster Pierre -Original Message- From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users- boun...@lists.osgeo.org] On Behalf Of Jason Mathis Sent: Tuesday, May 13, 2014 5:56 PM To: postgis-users@lists.osgeo.org Subject: [postgis-users] getting raster values for point from st_value Hi all, I have recently loaded a good amount of rasterized shape files into a db. I mainly used gdal_rasterize and raster2pgsql. The shape files were in NAD83 so I used the below commands to rasterized and load. gdal_rasterize -at -a desc -where desc='Very High' -burn 2 -tr 30 30 -ot byte -a_nodata -99 shape_file.shp raster.tiff raster2pgsql -d -s 5070 -t 100x100 -F -I -C -Y raster.tiff tbl_raster | psql -d raster_db I had two sets of shape files one set had about 10 burn values and another had about 50. The issue is getting these values out of the raster. I believe I should be using something from this page: http://postgis.net/docs/RT_ST_Value.html. I can visualize this file in QGIS and it looks good, the histogram shows values. I can also use the postgis functions, st_summarystats or st_histogram to see some data but I need to be able to query a given point and return the raster value. I have this: SELECT ST_Value(rast,(ST_Point( -111.750185, 34.886948))) FROM tbl_raster WHERE st_intersects(rast, (ST_Point(-111.750185, 34.886948))) Then i thought oh wait should I set the srid?: SELECT ST_Value(rast,(ST_SetSRID(ST_Point( -111.750185, 34.886948),5070))) FROM tbl_raster WHERE st_intersects(rast, (ST_SetSRID(ST_Point(-111.750185, 34.886948),5070))) Still I get nothing, then I thought I need to transform so I got this: SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_Point( -111.762866, 34.874309),5070),4326))) FROM tbl_raster WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_Point( - 111.762866, 34.874309),5070),4326))) But in the end nothing. I did need to add the SRID to the spatial_ref_sys table: INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 5070, 'EPSG', 5070, '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ', 'PROJCS[NAD83 / Conus Albers,GEOGCS[NAD83,DATUM[North_American_Datum_1983,SPHER OID[GRS 1980,6378137,298.257222101,AUTHORITY[EPSG,7019]],TOWGS84[0, 0,0,0,0,0,0],AUTHORITY[EPSG,6269]],PRIMEM[Greenwich,0,AUTHORIT Y[EPSG,8901]],UNIT[degree,0.0174532925199433,AUTHORITY[EPSG ,9122]],AUTHORITY[EPSG,4269]],PROJECTION[Albers_Conic_Equal_ Area],PARAMETER[standard_parallel_1,29.5],PARAMETER[standard_par allel_2,45.5],PARAMETER[latitude_of_center,23],PARAMETER[longitude _of_center,- 96],PARAMETER[false_easting,0],PARAMETER[false_northing,0],UNIT[ metre,1,AUTHORITY[EPSG,9001]],AXIS[X,EAST],AXIS[Y,NORTH],AUT HORITY[EPSG,5070]]’); But ultimately I have nothing I can’t seem to get any data out using a query with a point value. I am thinking it has to do with the srid? Since I can visually inspect it and see data. Anyone have any hot ideas as to why I am not seeing data? Thanks! This transmission contains confidential and privileged information intended solely for the party identified above. If you receive this message in error, you must not use it or convey it to others. Please destroy it immediately and contact the sender at (303) 386-3955 or by return e-mail to the sender. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] Postgis raster tutorial
Thank you Pierre, this will be of great help not just to me but to countless others as well, I hope it receives the widest possible circulation. -- View this message in context: http://postgis.17.x6.nabble.com/Postgis-raster-tutorial-tp5006124p5006287.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] getting raster values for point from st_value
Hi Pierre, Thanks so much for super quick response. You will have to forgive me I am fairly new to postgis. Creating the table is no problem but how would I check if it overlaps with my points? thanks! On May 13, 2014 at 4:08:40 PM, Pierre Racine (pierre.rac...@sbf.ulaval.ca) wrote: Create a new table with the extents of the raster tiles and check if it overlaps with your points: CREATE TABLE rastextets AS SELECT rast::geometry FROM tbl_raster Pierre -Original Message- From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users- boun...@lists.osgeo.org] On Behalf Of Jason Mathis Sent: Tuesday, May 13, 2014 5:56 PM To: postgis-users@lists.osgeo.org Subject: [postgis-users] getting raster values for point from st_value Hi all, I have recently loaded a good amount of rasterized shape files into a db. I mainly used gdal_rasterize and raster2pgsql. The shape files were in NAD83 so I used the below commands to rasterized and load. gdal_rasterize -at -a desc -where desc='Very High' -burn 2 -tr 30 30 -ot byte -a_nodata -99 shape_file.shp raster.tiff raster2pgsql -d -s 5070 -t 100x100 -F -I -C -Y raster.tiff tbl_raster | psql -d raster_db I had two sets of shape files one set had about 10 burn values and another had about 50. The issue is getting these values out of the raster. I believe I should be using something from this page: http://postgis.net/docs/RT_ST_Value.html. I can visualize this file in QGIS and it looks good, the histogram shows values. I can also use the postgis functions, st_summarystats or st_histogram to see some data but I need to be able to query a given point and return the raster value. I have this: SELECT ST_Value(rast,(ST_Point( -111.750185, 34.886948))) FROM tbl_raster WHERE st_intersects(rast, (ST_Point(-111.750185, 34.886948))) Then i thought oh wait should I set the srid?: SELECT ST_Value(rast,(ST_SetSRID(ST_Point( -111.750185, 34.886948),5070))) FROM tbl_raster WHERE st_intersects(rast, (ST_SetSRID(ST_Point(-111.750185, 34.886948),5070))) Still I get nothing, then I thought I need to transform so I got this: SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_Point( -111.762866, 34.874309),5070),4326))) FROM tbl_raster WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_Point( - 111.762866, 34.874309),5070),4326))) But in the end nothing. I did need to add the SRID to the spatial_ref_sys table: INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 5070, 'EPSG', 5070, '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ', 'PROJCS[NAD83 / Conus Albers,GEOGCS[NAD83,DATUM[North_American_Datum_1983,SPHER OID[GRS 1980,6378137,298.257222101,AUTHORITY[EPSG,7019]],TOWGS84[0, 0,0,0,0,0,0],AUTHORITY[EPSG,6269]],PRIMEM[Greenwich,0,AUTHORIT Y[EPSG,8901]],UNIT[degree,0.0174532925199433,AUTHORITY[EPSG ,9122]],AUTHORITY[EPSG,4269]],PROJECTION[Albers_Conic_Equal_ Area],PARAMETER[standard_parallel_1,29.5],PARAMETER[standard_par allel_2,45.5],PARAMETER[latitude_of_center,23],PARAMETER[longitude _of_center,- 96],PARAMETER[false_easting,0],PARAMETER[false_northing,0],UNIT[ metre,1,AUTHORITY[EPSG,9001]],AXIS[X,EAST],AXIS[Y,NORTH],AUT HORITY[EPSG,5070]]’); But ultimately I have nothing I can’t seem to get any data out using a query with a point value. I am thinking it has to do with the srid? Since I can visually inspect it and see data. Anyone have any hot ideas as to why I am not seeing data? Thanks! This transmission contains confidential and privileged information intended solely for the party identified above. If you receive this message in error, you must not use it or convey it to others. Please destroy it immediately and contact the sender at (303) 386-3955 or by return e-mail to the sender. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users -- This transmission contains confidential and privileged information intended solely for the party identified above. If you receive this message in error, you must not use it or convey it to others. Please destroy it immediately and contact the sender at (303) 386-3955 or by return e-mail to the sender. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
Re: [postgis-users] ST_MakeLine() woes
Remi, you replied before the edits I made to my last posting were accepted by the forum. It actually did work!, with very minor adjustments to your code. But I accept that my approach may not the best way to go about this problem, so I'll try your latest suggestion and let you know the results. Thank you for the new code. -- View this message in context: http://postgis.17.x6.nabble.com/ST-MakeLine-woes-tp5006266p5006289.html Sent from the PostGIS - User mailing list archive at Nabble.com. ___ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
