S-57 ENC coast lines, depth contours etc. imported incorrectly in OGR S57Reader

[Tom1]

S-57 reader in OGR incorrectly concatenates multi-part coast lines, depth contours, etc. found in S-57 ENC data. The problem is that features consisting of multiple polyline geometries are presented as one OGRLineString simply chaining the separate edges. This causes errors in ENC cell boundaries where e.g. coast line of an island cross the edge multiple times. (The coast line takes shortcuts to be continuous where it actually should be absent.) The fix involves creating a OGRMultiLineString geometry for such features that have multiple separate OGRLineString geometries embedded.

In the following source, there are added/removed tags for the changes required to released 1.5.0 source code.

Finally, here's the fixed source for S57Reader::AssembleLineGeometry?:

/************************************************************************/
/*                        AssembleLineGeometry()                        */
/************************************************************************/

void S57Reader::AssembleLineGeometry( DDFRecord * poFRecord,
                                     OGRFeature * poFeature )

{
   DDFField    *poFSPT;
   int         nEdgeCount;
   // OGRLineString *poLine = new OGRLineString();    // removed

/* -------------------------------------------------------------------- */
/*      Find the FSPT field.                                            */
/* -------------------------------------------------------------------- */
   poFSPT = poFRecord->FindField( "FSPT" );
   if( poFSPT == NULL )
       return;

   nEdgeCount = poFSPT->GetRepeatCount();

   OGRMultiLineString *poMultiLine =0; // added
   if(nEdgeCount>1) poMultiLine = new OGRMultiLineString(); // added

/* ==================================================================== */
/*      Loop collecting edges.                                          */
/* ==================================================================== */
   for( int iEdge = 0; iEdge < nEdgeCount; iEdge++ )
   {
       DDFRecord       *poSRecord;
       int             nRCID;

/* -------------------------------------------------------------------- */
/*      Find the spatial record for this edge.                          */
/* -------------------------------------------------------------------- */
       nRCID = ParseName( poFSPT, iEdge );

       poSRecord = oVE_Index.FindRecord( nRCID );
       if( poSRecord == NULL )
       {
           CPLError( CE_Warning, CPLE_AppDefined,
                     "Couldn't find spatial record %d.\n"
                     "Feature OBJL=%s, RCID=%d may have corrupt or"
                     "missing geometry.",
                     nRCID,
                     poFeature->GetDefnRef()->GetName(),
                     poFRecord->GetIntSubfield( "FRID", 0, "RCID", 0 ) );
           continue;
       }
  /* -------------------------------------------------------------------- */
/*      Establish the number of vertices, and whether we need to        */
/*      reverse or not.                                                 */
/* -------------------------------------------------------------------- */
       int             nVCount;
       int             nStart, nEnd, nInc;
       DDFField        *poSG2D = poSRecord->FindField( "SG2D" );
       DDFField        *poAR2D = poSRecord->FindField( "AR2D" );
       DDFSubfieldDefn *poXCOO=NULL, *poYCOO=NULL;

       if( poSG2D == NULL && poAR2D != NULL )
           poSG2D = poAR2D;

       if( poSG2D != NULL )
       {
           poXCOO = poSG2D->GetFieldDefn()->FindSubfieldDefn("XCOO");
           poYCOO = poSG2D->GetFieldDefn()->FindSubfieldDefn("YCOO");

           nVCount = poSG2D->GetRepeatCount();
       }
       else
           nVCount = 0;

       if( poFRecord->GetIntSubfield( "FSPT", 0, "ORNT", iEdge ) == 2 )
       {
           nStart = nVCount-1;
           nEnd = 0;
           nInc = -1;
       }
       else
       {
           nStart = 0;
           nEnd = nVCount-1;
           nInc = 1;
       }

/* -------------------------------------------------------------------- */
/*      Add the start node, if this is the first edge.                  */
/* -------------------------------------------------------------------- */
             OGRLineString *poLine = new OGRLineString(); // added

//        if( iEdge == 0 )        // removed
 //      {                        // removed
           int         nVC_RCID;
           double      dfX, dfY;
                     if( nInc == 1 )
               nVC_RCID = ParseName( poSRecord->FindField( "VRPT" ), 0 );
           else
               nVC_RCID = ParseName( poSRecord->FindField( "VRPT" ), 1 );

           if( FetchPoint( RCNM_VC, nVC_RCID, &dfX, &dfY ) )
               poLine->addPoint( dfX, dfY );
           else
               CPLError( CE_Warning, CPLE_AppDefined,
                         "Unable to fetch start node RCID%d.\n"
                         "Feature OBJL=%s, RCID=%d may have corrupt or"
                         " missing geometry.",
                         nVC_RCID,
                         poFeature->GetDefnRef()->GetName(),
                         poFRecord->GetIntSubfield( "FRID", 0, "RCID", 0 ) );
//        }    // removed
      /* -------------------------------------------------------------------- */
/*      Collect the vertices.                                           */
/* -------------------------------------------------------------------- */
       int             nVBase = poLine->getNumPoints();
             poLine->setNumPoints( nVCount+nVBase );

       for( int i = nStart; i != nEnd+nInc; i += nInc )
       {
           double      dfX, dfY;
           const char  *pachData;
           int         nBytesRemaining;

           pachData = poSG2D->GetSubfieldData(poXCOO,&nBytesRemaining,i);
                     dfX = poXCOO->ExtractIntData(pachData,nBytesRemaining,NULL)
                   / (double) nCOMF;
                     pachData = poSG2D->GetSubfieldData(poYCOO,&nBytesRemaining,i);
                     dfY = poXCOO->ExtractIntData(pachData,nBytesRemaining,NULL)
               / (double) nCOMF;
                         poLine->setPoint( nVBase++, dfX, dfY );
       }

/* -------------------------------------------------------------------- */
/*      Add the end node.                                               */
/* -------------------------------------------------------------------- */
       {
           int         nVC_RCID;
           double      dfX, dfY;
                     if( nInc == 1 )
               nVC_RCID = ParseName( poSRecord->FindField( "VRPT" ), 1 );
           else
               nVC_RCID = ParseName( poSRecord->FindField( "VRPT" ), 0 );

           if( FetchPoint( RCNM_VC, nVC_RCID, &dfX, &dfY ) )
               poLine->addPoint( dfX, dfY );
           else
               CPLError( CE_Warning, CPLE_AppDefined,
                         "Unable to fetch end node RCID=%d.\n"
                         "Feature OBJL=%s, RCID=%d may have corrupt or"
                         " missing geometry.",
                         nVC_RCID,
                         poFeature->GetDefnRef()->GetName(),
                         poFRecord->GetIntSubfield( "FRID", 0, "RCID", 0 ) );
       }
//    } // removed

       if( poLine->getNumPoints() >= 2 ){
//            poFeature->SetGeometryDirectly( poLine ); // removed
              if(poMultiLine) poMultiLine->addGeometryDirectly( poLine ); // added 
              else poFeature->SetGeometryDirectly( poLine );    // added
       }
       else
           delete poLine;
   } // added
     if(poMultiLine) poFeature->SetGeometryDirectly( poMultiLine );// added
}

[warmerdam]

Tom,

Would you have a redistributable file demonstrating the problem that this fixes?

[warmerdam]

S57 file where COALNE feature treated as multilinestring after patch.

[warmerdam]

Proposed Patch

[warmerdam]

Tom,

I have applied your change with a bit of reformatting in my "trunk" tree, and it seems to work. I found that the very small S-57 file in the autotest is affected by the change. In particular the COALNE feature geometry becomes a MULTILINESTRING instead of LINESTRING.

However, I notice that in the case of the COALNE feature, the geometry is actually a connected line string (now broken into two LINESTRINGs in a MULTILINESTRING). It seems to me that things should only be broken into a MULTILINESTRING when the lines can't be assembled into a contiguous linestring based on shared nodes.

I also notice that the layer type is not adjusted appropriately. If any line layer may return multilinestrings we will have to change their geometry types from wkbLineString to wkbUnknown. I can handle this aspect if needed.

So, currently I'm holding back on committing any change till we resolve whether contiguous linestrings can be composed into one.

[Tom1]

Replying to warmerdam:

Tom, Would you have a redistributable file demonstrating the problem that this fixes?

Unfortunately I do not have a redistributable file. (License of the ENC cell I have showing this phenomenon prohibits redistribution.)

Basically the geometry giving the trouble is trivial: For example, just a coast line of a small island having two capes projecting across the ENC cell boundary leaving a small part of the coast line (a small bay) inside the ENC cell will show the problem. Both parts should be separate (i.e. not connected) LINESTRINGs but still belonging to the same feature (the islands coast line).

[Tom1]

Replying to warmerdam:

Tom, I have applied your change with a bit of reformatting in my "trunk" tree, and it seems to work. I found that the very small S-57 file in the autotest is affected by the change. In particular the COALNE feature geometry becomes a MULTILINESTRING instead of LINESTRING. However, I notice that in the case of the COALNE feature, the geometry is actually a connected line string (now broken into two LINESTRINGs in a MULTILINESTRING). It seems to me that things should only be broken into a MULTILINESTRING when the lines can't be assembled into a contiguous linestring based on shared nodes. I also notice that the layer type is not adjusted appropriately. If any line layer may return multilinestrings we will have to change their geometry types from wkbLineString to wkbUnknown. I can handle this aspect if needed. So, currently I'm holding back on committing any change till we resolve whether contiguous linestrings can be composed into one.

Frank,

I must admit I'm not very familiar with the OGR internals and the requirements that comes from that side. I trust you are in position to evaluate and resolve the compatibility issues of the fix. For one, I understand from your reply that the layer geometry types need to be changed from wkbLineString to wkbUnknown. This is actually something with wider effect, since at least the depth contours have this same problem as the coast lines. Then again, I do not see any reason that other line features in S-57 files should be any different in this aspect. I guess any layer carrying LINESTRING geometries originated from S-57 should be able to carry MULTILINESTRING geometries as well.

I also understand you want to connect any chopped LINESTRINGs into one as far as possible. That sounds logical, but will obviously require a bit tinkering.

How should we proceed with making the fix that takes all the above factors into account? I'm afraid I'm not familiar enough with the code to make all the above fixes correctly. If you have the time to make them, I will test it with the ENCs I have.

Best regards,

Tom

[warmerdam]

Bogdan Grama has run into the same problem with a Romanian file. The file is public and at:

 http://download.osgeo.org/gdal/data/s57/bug2147_3R7D0889.000

Note that the COALNE feature in 1B5X02NE.000 should not be treated as a multilinestring, so it is not a good example.

I have reimplemented the patch to be careful to keep things as linestring even if there are multiple component lines as long as they connect up properly. I have also made a change to the layer type so that all object class line layers are now geometry type "unknown" since they might contain either linestring or multilinestring geometries. Urg.

Change committed in trunk (r19774) with some testing.

[warmerdam]

Change backported to 1.7 (r19975).