Make sure your layers align properly before importing them in the DB and don't forget to answer Guido's question...
> -----Original Message----- > From: Jason Mathis [mailto:[email protected]] > Sent: Wednesday, May 14, 2014 12:35 PM > To: Pierre Racine; PostGIS Users Discussion > Subject: RE: [postgis-users] getting raster values for point from st_value > > Its NAD83/Conus Albers > http://epsg.io/5070 > > I needed to reproject the shape file as well. It was originally in an Albers > projection but before handed over to me was projected into wgs84. The > raster and shape files did not line up properly until we reprojected the > shape file first and then did the conversion to raster. > > I had some issues with the reprojection as well and someone on the gdal > mailing list suggested using the epsg:5070 for the reprojection. So the > whole process really went > > 1. reproject shapfile > 2. create raster files > 3. Import files in the database > > I just reprojected one file using 4269 and it fits in the extents. Although > the > point is still hanging out there. > > I created the point like: > > SELECT ST_Transform(ST_SetSRID(ST_MakePoint(-111.750185, > 34.886948),5070),5070)::geometry(Point,5070) > > Which should be in arizona. > > > > On May 14, 2014 at 9:52:31 AM, Pierre Racine ([email protected]) > wrote: > > > What is this -s 5070 that you put in your raster2pgsql command? > Should not it be 4269 (NAD 83)? > > > -----Original Message----- > > From: Jason Mathis [mailto:[email protected]] > > Sent: Wednesday, May 14, 2014 11:47 AM > > To: Pierre Racine; PostGIS Users Discussion > > Subject: RE: [postgis-users] getting raster values for point from > st_value > > > > Well I was visualizing the tiff file not the postgis raster table. Its > seems > > somewhere along the lines the qgis raster plugin dropped off. > Although I am > > able to see the raster tables in the db manager and can drag/drop > the table, > > horribly slow but works. > > > > > > The raster extents or just the raster table do not overlay with the > point, > > which explains the problem. I still don’t understand what I did > wrong. > > Attaching a screen shot. > > > > thanks, > > > > > > > > On May 14, 2014 at 7:01:27 AM, Pierre Racine > ([email protected]) > > wrote: > > > > > > You said you could visualize this table in QGIS? Or were you > > speaking about a shapefile? You can visualize the rastextets table > and the > > points table the same way. > > > > Pierre > > > > > -----Original Message----- > > > From: Jason Mathis [mailto:[email protected]] > > > Sent: Tuesday, May 13, 2014 6:30 PM > > > To: PostGIS Users Discussion; Pierre Racine > > > Subject: 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 > > ([email protected]) > > > 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: [email protected] [mailto:postgis- > > > users- > > > > [email protected]] On Behalf Of Jason Mathis > > > > Sent: Tuesday, May 13, 2014 5:56 PM > > > > To: [email protected] > > > > 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 > > > [email protected] > > > 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. > > > > > > > > > > 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. > > > > > 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 [email protected] http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
