1

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

3 Answers 3

3

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.

1
  • 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, 2012 at 14:53
2

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.

4
  • 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, 2012 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/… Aug 1, 2012 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, 2012 at 14:56
  • 1
    Updated the answer above Aug 1, 2012 at 16:19
0

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.

2
  • Yes, I can do it in QGis, but I have to automate the process. It has to be done "in background".
    – exfizik
    Aug 1, 2012 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. Aug 1, 2012 at 11:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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