Monday, February 06, 2006

The Power of the Spatial to Attribute Join-Part 3

Sorry for the long delay between posts!

Today we are going to look at multi-dimensional subsurface mapping. First a little background. The Kansas Geological Survey has an extensive database of subsurface tops. A geologic top is the depth in feet from the surface to a particular geologic unit (system, group, formation, member). This data is collected from numerous scientists working over many years. When a well is drilled (mostly for petroleum) geophysical logs are run in the hole to measure different properties of the rocks. These electric logs are scanned, studied, and geologic units are picked from the logs (Chase group, Permian system, Fort Scott limestone formation). Read this primer for more information.

So what does this mean for GIS and SQL? With this database we have access to over 1,700,000 called tops spread across ~150,000 wells in the state of Kansas. Each well can/will have more than one called formation. We have more than one record per geographic location--in this case the records represent mulitple entries (formations) in the z domain.

We can visualize some of this dataset in ArcScene (or any other 3D package).
Most analysis and visualization occurs when users map one unit of measure (depth to unit, thickness between units, or something else entirely) per geographic location.

The query to select an individual unit is pretty easy.

SQL> desc strat_well_top_1JUN2005
Name Null? Type
--------------------------------- -------- ------------------------
FORMATION_NAME NOT NULL VARCHAR2(60)
WELL_HEADER_KID NOT NULL NUMBER(10)
WELL_TOP NUMBER(8,2)
FIELD_KID NUMBER(10)
LATITUDE NUMBER(11,6)
LONGITUDE NUMBER(11,6)
GROUND_ELEVATION VARCHAR2(40)
GROUND_ELEVATION_SOURCE VARCHAR2(12)
WELL_TOP_SEALEVEL_NON_NED VARCHAR2(40)
WELL_TOP_SEALEVEL NUMBER

1 select well_header_kid, latitude, longitude,
2 well_top, well_top_sealevel
3 from strat_well_top_1JUN2005
4* where FORMATION_NAME = 'Permian System'
SQL> /


WELL_HEADER_KID LATITUDE LONGITUDE WELL_TOP WELL_TOP_SEALEVEL
--------------- ---------- ---------- ---------- -----------------
1025687130 37.58751 -101.36115 680 -2381
1028094976 37.61659 -100.73941 834 -2095
1006155619 37.39194 -101.67532 608 -2673
1006514758 39.99071 -97.14744 358 -1158
1025773538 37.20899 -101.90656 1395 -2120
1024492073 37.00136 -100.56991 3550 1170
1026599879 37.91949 -101.48087 660 -2529
1006054521 37.50601 -99.70401 387 -2066
1006054213 37.53783 -99.8804 417 -2098
1006054635 37.48215 -100.02145 412 -2179
1002915548 37.55853 -100.38357 775 -2008


This can be done easily outside of database. With ArcMap you can do a definition query to limit the tabluar dataset to only the unit in question.

With ArcIMS you can perform an ArcXML spatial query to only select the unit in question.
<LAYER type="featureclass" name="Depth" visible="true">
<DATASET fromlayer="1" />
<SPATIALQUERY searchorder="attributefirst"
where="PLSS.WELL_HEADERS_29AUG2003.KID = PLSS.STRAT_WELL_TOP_28AUG2003.WELL_HEADER_KID
and PLSS.STRAT_WELL_TOP_28AUG2003.FORMATION_NAME = 'Permian System'
and PLSS.STRAT_WELL_TOP_28AUG2003.WELL_TOP <> 9999
and PLSS.STRAT_WELL_TOP_28AUG2003.GROUND_ELEVATION > 0"
jointables="PLSS.STRAT_WELL_TOP_28AUG2003" />
<VALUEMAPRENDERER lookupfield="PLSS.STRAT_WELL_TOP_28AUG2003.WELL_TOP">
<RANGE lower="15" upper="326" label="15 TO 326">
<SIMPLEMARKERSYMBOL color="56,168,0" width="6" />
</RANGE>
<RANGE lower="326" upper="637" label="326 TO 637">
<SIMPLEMARKERSYMBOL color="90,186,0" width="6" />
</RANGE>
<RANGE lower="637" upper="949" label="637 TO 949">
<SIMPLEMARKERSYMBOL color="131,207,0" width="6" />
</RANGE>
<RANGE lower="949" upper="1260" label="949 TO 1260">
<SIMPLEMARKERSYMBOL color="176,224,0" width="6" />
</RANGE>
<RANGE lower="1260" upper="1571" label="1260 TO 1571">
<SIMPLEMARKERSYMBOL color="228,245,0" width="6" />
</RANGE>
<RANGE lower="1571" upper="1882" label="1571 TO 1882">
<SIMPLEMARKERSYMBOL color="255,225,0" width="6" />
</RANGE>
<RANGE lower="1882" upper="2193" label="1882 TO 2193">
<SIMPLEMARKERSYMBOL color="255,170,0" width="6" />
</RANGE>
<RANGE lower="2193" upper="2504" label="2193 TO 2504">
<SIMPLEMARKERSYMBOL color="255,115,0" width="6" />
</RANGE>
<RANGE lower="2504" upper="2816" label="2504 TO 2816">
<SIMPLEMARKERSYMBOL color="255,55,0" width="6" />
</RANGE>
<RANGE lower="2816" upper="3550" label="2816 TO 3550">
<SIMPLEMARKERSYMBOL color="255,0,0" width="6" />
</RANGE>
</VALUEMAPRENDERER>
</LAYER>

This looks like....


Lets get a little more complicated. Say you want to compute the thickness between two different formations. One way to do this is to select the two formations of interest, group the results, and select the difference between the min depth value and the max depth value.


select WELL_HEADER_KID, round(max(WELL_TOP) - min(WELL_TOP),1) STRAT_THICKNESS
from strat_well_top_1JUN2005
where (formation_name = 'Chase Group'
or formation_name = 'Council Grove Group')
and WELL_TOP is not null
and WELL_TOP <> 9999
and LONGITUDE >= -103.57371
and LONGITUDE <= -98.352965
and LATITUDE >= 34.9908216
and LATITUDE <= 40.2115692
group by WELL_HEADER_KID
having count(distinct formation_name) = 2
and max(WELL_TOP) - min(WELL_TOP) > 0
and max(WELL_TOP) - min(WELL_TOP) <> max(WELL_TOP)


By selecting only the Chase Group and the Council Grove Group we lower our tops table to records that have either the Chase or The Council Grove. Once we group by the well_header_kid and select records that have rows in both groups we can do the math to compute the thickness.
group by WELL_HEADER_KID 
having count(distinct formation_name) = 2

This drops records of wells that have one group or the other, but not both.

Now we can subtract the "min(well_top)" depth from the "max(well_top)" depth to compute the thickness between the two groups for a particluar well.

Slap that query to a view and join that in ArcIMS and you have nice point map of the thickness between two geologic groups.


Give it a try at: http://hercules.kgs.ku.edu/kgs/oilgas/strat_welltops/top_viewer.cfm

Or we can use ArcGIS server to generate a raster dataset off that that view (example application coming soon).

What do you think?

Cheers,

Jeremy

Comments:
Hey Jeremy,
I'm wondering if you know how to add an SQL REMARK to an ArcMAP definition query. I would like to use something like one of the following:
REMARK
--
/* */
None of these seem to allow to the definiton guery to execute. Any ideas? Please let me know when you get a chance.

Regards,
thegis.blogspot.com
 
Good question. You can do C++ style comments in the definition query.

/* test */
id > 2

It worked for me at least. Just remember to close the wrap the text with a beginning /* and an ending */.

This is not a standard SQL remark because ArcMap definition queries are only for filtering data. It works for most data formats that ArcGIS supports, but you must use nomenclature that matches the data format that you are querying. Wildcard for sde: "%"
Wildcard for personal geodatabase: "*"

Cheers,

Jeremy
 
Thanks for the quick response!
Unfortunately I had tried using
/* Comments */
already and it didn't work. It doesn't complain but no results are displayed. If I remove the /* Comments */ then the results are shown...
It would be nice to be able to do this since sometimes I go back to layerfiles or mxds and wonder what the heck I was querying for?

Cheers,
daniel
 
Hey Daniel,

That is interesting. I am successful if I do the following definition query on a table from Oracle.

/* comment */
FIELD_KID = 1000146432

I am using ArcGIS 9.1 and executing this query on an event theme from a table that lives in Oracle.

If I try to do the following on a shapefile I get an error.

/* comment for shapefile */
"CNTY_FIPS" = '007'

It looks like that you are only able to do this on certain types of data formats. What data type are you trying this on?

The real issue it seems is that you would like to be able save metadata with your layer file. Currently you would have to plug that metadata into the "Description" text field under the General tab of the layer properties.

I think that ESRI is working on a more encompassing geodatabase for metadata, layer files, models, etc. Probably post 9.2 though.

Cheers,

Jeremy
 
Thanks again Jeremy!

Yeah I've tried it unsuccessfully with shape, cov, personal gdb and SDE data types.

We are using ArcGIS 9.1 - hopefully 9.2 will provide more layerfile documenting capabilities. I liked your idea of adding a description to the layerfile though.

Cheers,
Daniel
 
Post a Comment



<< Home

This page is powered by Blogger. Isn't yours?