Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a bunch of ESRI shapefiles that contain Land Cover data across Canada (provided by Geobase). I need to query the files and obtain the value of a certain field (called COVTYPE) in a given location specified either by lat/lon or UTM coordinates. I can open a file in Quantum GIS and find the needed values manually, but this process needs to be done automatically for a large number of values so I'm using C++ and GDAL.

I'm very new to the whole GIS area and I'm just starting to learn GDAL. So far I've learned how to loop over layers, features, fields etc, but I can't figure out how to relate coordinates to features. I went through the official GDAL/OGR tutorials, but I'm still stuck. Any pointers would be appreciated!

Update: Just to clarify. My input data would be coordinates (lat/lon or UTM). So the way I imagine this should work is:

  1. I find which feature the given point is located within
  2. read the field value in that feature

I know how to do step 2, but I'm stuck at step 1.

Final solution Technically, the update to Ragi's answer comes very close to what I need to do, but I'll post here a quick summary of how I ended up using it. My way is easier for constructing a point geometry object.

So leaving all the boilerplate code out, the essential part is:

OGRPoint* point = new OGRPoint(some_longitude, some_latitude);
layer->SetSpatialFilter(point); //getting only the feature intercepting the point
OGRFeature* feature = layer->GetNextFeature();
int value = feature->GetFieldAsInteger(0); //getting the value - DONE!
//in my case it's 1st field, hence the 0th index
share|improve this question

3 Answers 3

up vote 2 down vote accepted

If you look at the OGR API Tutorial you will find this piece of code which does something very close to what you want:

#include "ogrsf_frmts.h"

int main()
{
    OGRRegisterAll();

    OGRDataSource       *poDS;

    poDS = OGRSFDriverRegistrar::Open( "point.shp", FALSE );
    if( poDS == NULL )
    {
        printf( "Open failed.\n" );
        exit( 1 );
    }

    OGRLayer  *poLayer;

    poLayer = poDS->GetLayerByName( "point" );

    OGRFeature *poFeature;

    poLayer->ResetReading();
    while( (poFeature = poLayer->GetNextFeature()) != NULL )
    {
        OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
        int iField;

        for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
        {
            OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );

            if( poFieldDefn->GetType() == OFTInteger )
                printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
            else if( poFieldDefn->GetType() == OFTReal )
                printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
            else if( poFieldDefn->GetType() == OFTString )
                printf( "%s,", poFeature->GetFieldAsString(iField) );
            else
                printf( "%s,", poFeature->GetFieldAsString(iField) );
        }

        OGRGeometry *poGeometry;

        poGeometry = poFeature->GetGeometryRef();
        if( poGeometry != NULL 
            && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
        {
            OGRPoint *poPoint = (OGRPoint *) poGeometry;

            printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() );
        }
        else
        {
            printf( "no point geometry\n" );
        }       
        OGRFeature::DestroyFeature( poFeature );
    }

    OGRDataSource::DestroyDataSource( poDS );
}

The relevant lines are

poFeature->GetFieldAsInteger( iField )
poFeature->GetFieldAsDouble(iField)
poFeature->GetFieldAsString(iField)

Update: Based on the comments below, it sounds like what you need is a spatial query. To do a spatial query, all you need is to create a geometry object and set it before the loop. The easiest way I can think of creating your query geometry is to simply use one of the methods in OGRGeometryFactory (like OGRGeometryFactory::ConstructFromWKT) to the a OGRGeometry object, then you pass that object to OGRLayer::SetSpatialFilter.

There are plenty of examples in the OGR test suite. Nevertheless, the ones that are easier to understand (i.e. that don't have much of test harness cruft) are the python tests. I would start with those. The API is very symmetrical to the C++ one, so translating from one to the other should be very easy.

share|improve this answer
    
Yes, I've gotten that far - I know how to get the required field values for a certain feature in a certain layer (there is actually just 1 layer in these files). But how do I know which feature I need to read? –  exfizik Aug 1 '12 at 2:42
1  
Are you asking how to figure out which point intersects a particular polygon? You can do that with a Spatial Query. Before you start looping, you set a spatial filter gdal.org/ogr/… –  Ragi Yaser Burhum Aug 1 '12 at 4:17
    
Yes, I think it sounds like what I need to do. I'll give it a try and see where it leads me. Do you happen to have any code examples for spatial queries? –  exfizik Aug 1 '12 at 14:56
1  
Updated the answer above –  Ragi Yaser Burhum Aug 1 '12 at 16:19
    
Thanks for the pointer! –  exfizik Aug 3 '12 at 20:09

To find which feature the coordinate is near, you first need to build an R-Tree index of bounding boxes or envelope of each feature. A popular library for this is libspatialindex.

Secondly, you would then need to know for each of the matched features from your R-tree, which ones match. GDAL/OGR does have some operations based on GEOS to see if the geometry "contains" the point of interest, and you can extract the field info. See the OGRGeometry Class Reference for Contains, Touches, Within, etc. that can perform the geometry relation operator.


Personally I wouldn't do any of the above, because it would take me too long to figure out the coding, and I know there are faster ways. (I would loading all shapefiles in a PostGIS database, then query the locations using SQL.) However, it is understandable if this GIS functionality needs to be embedded in an existing project without adding more complicated dependencies, such as a database server.

share|improve this answer
    
Thanks for the ideas. I have to do it in C++ only - adding a database to this is not an option at this point. It probably would be easier if I had all the data in a DB, but right now I have a bunch of separate files and it's likely to stay that way in the foreseeable future. –  exfizik Aug 1 '12 at 14:53

I'm not sure if/how you can do this in GDAL (or OGR technically, for vectors). However, you most certainly can do it easily in QGIS. First, create a point shapefile from your locations (if it's not already). Then you want to do a "Join Attributes By Location" available under the Vector > Data Management Tools Menu.

share|improve this answer
    
Yes, I can do it in QGis, but I have to automate the process. It has to be done "in background". –  exfizik Aug 1 '12 at 2:45
    
Ah, sorry. I didn't realize it had to be a 'background' job. In that case, I second Mike Toews PostGIS suggestion. –  Darren Cope Aug 1 '12 at 11:23

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.