GDAL/OGR documentation problem (S-57)

General Usage
1

Hello,

I have a Python-program to analyse data and I am rewriting it in Julia. The Python version needs over 14 hours, Julia will probably end in under 2 hours.

Now my problem. The result needs to be plotted on a maritime map, S-57 format. In Python I wrote a small library for the maps using GDAL/OGR. Julia has a GDAL package, but the documentation is far too limited for me: I cannot get anything done. Even opening a S-57 file fails in Julia.

Can someone point me to documentation that fits my need, or even help me on my way? I attach a bare minimum file in Python that has my most needed elements.

Thank you for any help.

from osgeo import ogr

def load_objects(infile, name):
    loc=[]
    dataset = ogr.Open( infile )
    try:
        objects=dataset.GetLayerByName( name )
    except:
        return loc
    if objects is not None:
        feat = objects.GetNextFeature()
        while feat is not None:
            geom=feat.GetGeometryRef()
            if geom.GetGeometryType() == ogr.wkbLineString:
                pnt_count = geom.GetPointCount()
                for pnt in range(pnt_count):
                    lat,lon= geom.GetY(pnt), geom.GetX(pnt)
                    loc.append((lat,lon))
    return loc

loc = load_objects(<some S-57 datafile>, "NAVLNE")
print(loc)
2

I made a step forward. GDAL may not work, ArchGDAL works:

dataset = ArchGDAL.read(infile)
objects = ArchGDAL.getlayer(dataset, name)

Now I need to find and extract an array of lat/lon’s.
So far I could find the methods needed:

API Reference · ArchGDAL.jl

Any suggestions?

3

OK, very few seem to have experience with my problem.
This is what happens:

julia> dataset = ArchGDAL.read(infile)
GDAL Dataset (Driver: S57/IHO S-57 (ENC))
File(s):
/home/ removed .000

Number of feature layers: 37
Layer 0: DSID (wkbNone)
Layer 1: BOYCAR (wkbPoint)
Layer 2: BOYLAT (wkbPoint)
Layer 3: BOYSAW (wkbPoint)
Layer 4: BOYSPP (wkbPoint)
Remaining layers:
CBLSUB, CTNARE, DWRTPT, DEPARE, DEPCNT,
EXEZNE, FOGSIG, LIGHTS, MAGVAR, MIPARE,
OBSTRN, OFSPLF, PIPSOL, PRCARE, RTPBCN,
RDOSTA, RESARE, SEAARE, SBDARE, SOUNDG,
TS_PAD, TOPMAR, TSSBND, TSSCRS, TSSLPT,
TSEZNE, WRECKS, M_COVR, M_NPUB, M_NSYS,
M_QUAL, C_AGGR,

julia> name = “DEPCNT”
“DEPCNT”

julia> objects = ArchGDAL.getlayer(dataset, name)
Layer: DEPCNT
Geometry 0 (): [wkbUnknown], LINESTRING (3.294656…), …
Field 0 (RCID): [OFTInteger], 274, 275, 276, 277, 278, 279, 280, 281, …
Field 1 (PRIM): [OFTInteger], 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
Field 2 (GRUP): [OFTInteger], 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
Field 3 (OBJL): [OFTInteger], 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, …
Field 4 (RVER): [OFTInteger], 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
…
Number of Fields: 24

julia> feat = ArchGDAL.unsafe_getfeature(objects, 1)
NULL Feature

How do I get to an array of lat/lon’s? (seems to be in Geometry 0) Later on I will need the depth-value too.
Thank you