Python OGR CreateGeometryFromGML does not handle unit for radius in CURVEPOLYGON/CIRCULARSTRING geometry

[dearith]

Expected behavior and actual behavior.

Python OGR CreateGeometryFromGML does not handle unit for radius in CURVEPOLYGON/CIRCULARSTRING geometry

Steps to reproduce the problem.

gml = '<gml:CircleByCenterPoint><gml:pos>49 2</gml:pos><gml:radius uom="[nmi_i]">2000</gml:radius></gml:CircleByCenterPoint>'
geom1 = ogr.CreateGeometryFromGML(gml)
print(geom1)

CIRCULARSTRING (-1951 2,2049 2,-1951 2)

gml = '<gml:CircleByCenterPoint><gml:pos>49 2</gml:pos><gml:radius uom="m">2000</gml:radius></gml:CircleByCenterPoint>'
geom1 = ogr.CreateGeometryFromGML(gml)
print(geom1)

CIRCULARSTRING (-1951 2,2049 2,-1951 2)

gml = '<gml:CircleByCenterPoint><gml:pos>49 2</gml:pos><gml:radius>2000</gml:radius></gml:CircleByCenterPoint>'
geom1 = ogr.CreateGeometryFromGML(gml)
print(geom1)

CIRCULARSTRING (-1951 2,2049 2,-1951 2)

Whatever the unit, the geometry is always the same.

Notice : if I try to display my GML file with QGIS, the uom of the radius element is handled by QGIS;

Operating system

CentOS 8, Windows 10

GDAL version and provenance

3.0.2 and 3.3.3 : CentOS 8 from http://download.osgeo.org/gdal/
3.0.2 and 3.3.3 : Windows10 from https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal

Python version : 3.7

[nyalldawson]

Notice : if I try to display my GML file with QGIS, the uom of the radius element is handled by QGIS;

This is odd -- QGIS uses OGR to read GML files, so you should see the same results outside of QGIS (in other words, QGIS doesn't have any specif handling of the uom attribute).

Do you have a sample GML file you were using to test that I can investigate?

[dearith]

@nyalldawson : thanks for you answer. Here is a GML file I'm using, in the attached zip file
LY_TEST_XPATH.zip

[jratike80]

Does that GML tell what is the CRS of CenterPoint? Just wondering that if GDAL thinks that the coordinates are in same units than uom the results would make sense. I believe that intention is that CenterPoint is in EPSG:4326, right?

[dearith]

Yes. The CRS is in EPSG:4326.

[jratike80]

I apologize that I do not master GML but could you try to edit the GML and tell the srsName of Pos here:

<gml:CircleByCenterPoint interpolation="circularArcCenterPointWithRadius" numArc="1">
           <gml:pos>-10.8 121.9</gml:pos>
           <gml:radius uom="[nmi_i]">140.0</gml:radius>
</gml:CircleByCenterPoint>

But probably it is OK because this appears earlier:
<aixm:Surface srsName="http://www.opengis.net/def/crs/EPSG/0/4326" srsDimension="2" axisLabels="Lat Long" gml:id="uuid.0f8f949f-9e6a-46c6-87ae-1a47b7aabd0a">

Corresponding unit test in the autotest suite is here:
https://github.com/OSGeo/gdal/blob/master/autotest/ogr/ogr_gml_geom.py#L1746

I believe that your test is indeed not good but my guess that srsName could be told for Pos was wrong. Or perhaps it could work as well, GML tend to have many ways for everything. But try to give the srsName directly for CircleByCenterPoint.
gml = '<gml:CircleByCenterPoint srsName="urn:ogc:def:crs:EPSG::4326"><gml:pos>49 2</gml:pos><gml:radius uom="[nmi_i]">2000</gml:radius></gml:CircleByCenterPoint>'
Perhaps the issue is if the srsName is correctly transferred to CircleByCenterPoint from the outer element aixm:Surface.

[dearith]

I have also try to set srsName, but setting the srsName doesn't change anything in my case, and then a warning is displayed:

gml = '<gml:CircleByCenterPoint srsName="urn:ogc:def:crs:EPSG::4326"><gml:pos>49 2</gml:pos><gml:radius uom="[nmi_i]">2000</gml:radius></gml:CircleByCenterPoint>'
geom1 = ogr.CreateGeometryFromGML(gml)
print(geom1)

ERROR 1: PROJ: proj_create: no database context specified
CIRCULARSTRING (-1951 2,2049 2,-1951 2)

[jratike80]

Did you do from osgeo import osr?

[dearith]

No, I don't import osr, because I don't use it.

from osgeo import ogr
gml = '<gml:CircleByCenterPoint srsName="urn:ogc:def:crs:EPSG::4326"><gml:pos>49 2</gml:pos><gml:radius uom="[nmi_i]">2000</gml:radius></gml:CircleByCenterPoint>'
geom1 = ogr.CreateGeometryFromGML(gml)
print(geom1)

Same result if I import osr.

from osgeo import ogr, osr
gml = '<gml:CircleByCenterPoint srsName="urn:ogc:def:crs:EPSG::4326"><gml:pos>49 2</gml:pos><gml:radius uom="[nmi_i]">2000</gml:radius></gml:CircleByCenterPoint>'
geom1 = ogr.CreateGeometryFromGML(gml)
print(geom1)

[jratike80]

Works for me

>>> print(geom1)
LINESTRING (82.3333333333333 2.0,81.9649708951405 17.9162820740288,80.9503899734996 31.0922331089589,79.4774053714158 40.7263319998154,77.7134143645019 47.3788749910705,75.770646747763 51.8719230075544,73.717960773297 54.862804779185,71.5975258680789 56.8045206362836,69.4360541467574 57.9982...

[dearith]

How do you do to have a CIRCULARSTRING as input and the have a LINESTRING as output ?
Notice, I don't use the gdal/ogr package from QGIS. I use
3.0.2 and 3.3.3 : CentOS 8 from http://download.osgeo.org/gdal/
3.0.2 and 3.3.3 : Windows10 from https://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal

[jratike80]

I do not know, I am just another GDAL user. Probably that is an intentional thing. Something about how arcs are handled is documented in https://gdal.org/drivers/vector/gml.html. For further questions you should use gdal-dev mailing list.

I somehow believe that the reason for creating a linear geometry in this case with ogr.CreateGeometryFromGML(gml) is in the combination of having the center point in EPSG:4326 and radius in linear units like nautical miles. Because 2000 nm corresponds to different number of degrees in north-south and west-east directions the result, when looked as EPSG:4326 coordinates, does not really look like a circle (I believe that the image below should be rotated because my geometry viewer is not aware of lat-long coordinate order but tilt your head and look at the shape).

[jratike80]

This one is related #3118. That fix was targeted into GDAL 3.2.1 so your GDAL 3.3.3 should support also nautical miles.

[dearith]

Yes, you are right. If srsName is set and using GDAL 3.3.3, it works. Thanks for your help.

[rouault]

If srsName is set and using GDAL 3.3.3, it works.
confirmed:

from osgeo import ogr
gml = '<gml:CircleByCenterPoint srsName="urn:ogc:def:crs:EPSG::4326"><gml:pos>49 2</gml:pos><gml:radius uom="[nmi_i]">2000</gml:radius></gml:CircleByCenterPoint>'
geom1 = ogr.CreateGeometryFromGML(gml)
print(geom1)
LINESTRING (82.3333333333333 2.0,81.9649708951405 17.9162820740288,80.9503899734996 31.0922331089589,79.4774053714158 40.7263319998154,77.7134143645019 47.3788749910705,75.770646747763 51.8719230075544,73.717960773297 54.862804779185,71.5975258680789 56.8045206362836,69.4360541467574 57.9982006490851,67.251283798058 58.6446279627804,65.0556159618402 58.8799722918398,62.8581818816355 58.7985588760807,60.6660509148471 58.4671384061892,58.484959403598 57.9339009526058,56.319764743939 57.23426907228,54.1747369589143 56.3946985802436,52.0537512120265 55.4352260434098,49.9604181222299 54.3712124513868,47.8981738886035 53.2145610052826,45.8703436964326 51.9745841841072,43.8801868351735 50.6586326318515,41.9309289037329 49.2725595553713,40.0257845818302 47.8210697576958,38.1679732423128 46.3079866069935,36.3607289000962 44.7364598599306,34.6073054811915 43.109130324446,32.9109780555496 41.4282626333621,31.2750404526033 39.6958541386365,29.7027995329023 37.9137256383827,28.1975663009263 36.083598000658,26.7626439990047 34.2071575446913,25.4013133111614 32.2861121474305,24.1168148224037 30.3222393728771,22.9123289187483 28.3174274161271,21.7909533719019 26.2737092764657,20.7556789256889 24.1932903011585,19.8093632841697 22.0785690597146,18.9547039882947 19.9321514087449,18.1942107524631 17.7568575841583,17.5301779074635 15.5557222047693,16.9646576546731 13.3319871824484,16.4994348711205 11.0890876991186,16.1360042100484 8.83063161719208,15.8755502126393 6.56037292111702,15.7189310815552 4.28218002455675,15.6666666666667 2.0,15.7189310815552 -0.282180024556753,15.8755502126393 -2.56037292111702,16.1360042100484 -4.83063161719208,16.4994348711205 -7.08908769911865,16.9646576546731 -9.33198718244836,17.5301779074635 -11.5557222047693,18.1942107524631 -13.7568575841583,18.9547039882947 -15.9321514087449,19.8093632841697 -18.0785690597146,20.7556789256889 -20.1932903011585,21.7909533719019 -22.2737092764657,22.9123289187483 -24.3174274161271,24.1168148224037 -26.3222393728771,25.4013133111614 -28.2861121474305,26.7626439990047 -30.2071575446913,28.1975663009263 -32.083598000658,29.7027995329023 -33.9137256383827,31.2750404526033 -35.6958541386365,32.9109780555496 -37.4282626333621,34.6073054811915 -39.109130324446,36.3607289000962 -40.7364598599306,38.1679732423128 -42.3079866069935,40.0257845818302 -43.8210697576958,41.9309289037329 -45.2725595553713,43.8801868351735 -46.6586326318515,45.8703436964326 -47.9745841841072,47.8981738886035 -49.2145610052826,49.9604181222299 -50.3712124513868,52.0537512120265 -51.4352260434098,54.1747369589143 -52.3946985802436,56.319764743939 -53.23426907228,58.484959403598 -53.9339009526058,60.666050914847 -54.4671384061892,62.8581818816355 -54.7985588760807,65.0556159618402 -54.8799722918398,67.251283798058 -54.6446279627804,69.4360541467574 -53.9982006490851,71.5975258680789 -52.8045206362836,73.717960773297 -50.862804779185,75.770646747763 -47.8719230075544,77.7134143645019 -43.3788749910705,79.4774053714158 -36.7263319998154,80.9503899734995 -27.0922331089592,81.9649708951405 -13.9162820740288,82.3333333333333 2.0)

[vlapra]

It works if the srsName is set directly on the CircleByCenterPoint element (returns a POLYGON). It also works if the srsName is defined on some parent GML elements, but not on all.

However, if the srsName is only defined on the root element (Surface), it doesn't work (returns a CURVEPOLYGON).

Is there any workaround? When processing large files, modifying the XML to insert the correct srsName for child elements is inconvenient.

PS: It is strange that when I replace nmi_i with m, km, or nm, it doesn't work (regardless of the presence of srsName) — it returns a CURVEPOLYGON. According to the comments above and #3118, it should work.

<aixm:Surface gml:id="ID_240" srsName="urn:ogc:def:crs:EPSG::4326">
  <gml:patches>
    <gml:PolygonPatch>
      <gml:exterior>
        <gml:Ring>
          <gml:curveMember>
            <gml:Curve gml:id="ID_241">
        	<gml:segments>
        	  <gml:CircleByCenterPoint numArc="1" srsName="urn:ogc:def:crs:EPSG::4326">
        	    <gml:posList>52.36666666666667 -22.1</gml:posList>
        	    <gml:radius uom="[nmi_i]">15.0</gml:radius>
        	  </gml:CircleByCenterPoint>
                </gml:segments>
             </gml:Curve>
          </gml:curveMember>
        </gml:Ring>
     </gml:exterior>
   </gml:PolygonPatch>
 </gml:patches>
 <aixm:horizontalAccuracy uom="KM">2.0</aixm:horizontalAccuracy>
</aixm:Surface>

GDAL 3.11.0 (build from master branch)