Checking Type Locations
Submitted by dylan on Mon, 2009-04-20 22:18.
Just Checking
-- NAD27 to NAD83
echo 119d7\'4\"W 36d23\'13\"N | cs2cs +proj=latlong +datum=NAD27 +to +proj=latlong +datum=NAD83 -f "%.6f"
-119.118718 36.386894 0.000037
-- the relevent SRIDs
srid | proj4text
------+----------------------------------------------------
4267 | +proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs
4269 | +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs
4326 | +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
-- NAD27 type location point
\SET pt Transform(SetSRID(ST_MakePoint(-119.1177777, 36.38694444), 4267), 9001)
-- find intersecting map unit, and the NAD83 coords
SELECT AsText(Transform(:pt, 4269)) AS NAD83_pt, mapunit_poly.mukey, muname
FROM mapunit_poly
JOIN mapunit
ON mapunit_poly.mukey = mapunit.mukey
AND st_intersects(wkb_geometry, :pt) ;
nad83_pt | mukey | muname
-------------------------------------------+--------+--------------------------------------------------------
POINT(-119.118718055326 36.3868941884307) | 467194 | Grangeville sandy loam, drained, 0 TO 2 percent slopes
-- 250 m buffer around NAD27 type location coords
\SET box st_buffer(Transform(SetSRID(ST_MakePoint(-119.1177777, 36.38694444), 4267), 9001), 250)
-- select overlapping map units and their associated percent overlap
SELECT mapunit_poly.mukey, muname,
round((ST_Area(ST_Intersection(wkb_geometry, :box)) / ST_Area(:box)) * 100) AS area_pct
FROM mapunit_poly
JOIN mapunit
ON mapunit_poly.mukey = mapunit.mukey
AND st_intersects(wkb_geometry, :box)
ORDER BY area_pct DESC ;
mukey | muname | area_pct
--------+--------------------------------------------------------+----------
467194 | Grangeville sandy loam, drained, 0 TO 2 percent slopes | 70
463596 | Grangeville silt loam, drained | 12
467206 | Riverwash | 11
467210 | Tujunga loamy sand, 0 TO 2 percent slopes | 7
|
