1

我有一个包含多个多边形的 shapefile，这是我的文件的简短版本（简单功能）：

mygeom<-structure(list(itemID = c("a1234", "a27"), status = c(1, 0), 
    Schema1 = c("A", "A"), Schema2 = c("B", NA), Date = c("2020-03-01", 
    "2020-03-01"), geometry = structure(list(structure(list(structure(c(-0.649062881756964, 
    -0.649112159426678, -0.649217037157584, -0.649329802261076, 
    -0.649385577124541, -0.64940943118353, -0.649471971919662, 
    -0.649418938777326, -0.649353483653088, -0.64927763496891, 
    -0.64916144166702, -0.649060881454067, -0.648865049171286, 
    -0.648768067415234, -0.648597397154666, -0.648913174239608, 
    -0.648976607280396, -0.649027256609277, -0.649058695235838, 
    -0.649062881756964, 41.4921932026922, 41.4922493593042, 41.4923567934804, 
    41.4924515944872, 41.492509205028, 41.4925254105761, 41.4925678967614, 
    41.4926030387555, 41.4926842652395, 41.4927413468292, 41.4927247423209, 
    41.4926502665709, 41.4924908972756, 41.4924318736095, 41.4923280025819, 
    41.4920576807899, 41.4921104610734, 41.4921526044294, 41.4921884317747, 
    41.4921932026922), .Dim = c(20L, 2L))), class = c("XY", "POLYGON", 
    "sfg")), structure(list(structure(c(-0.886732800166245, -0.886458452915082, 
    -0.886077569154315, -0.885592339912652, -0.885526044922013, 
    -0.885468424362897, -0.885411911891854, -0.885182818538855, 
    -0.8850352039955, -0.885175736668516, -0.885147824091142, 
    -0.885184202391753, -0.885184872943689, -0.885184956762793, 
    -0.885186297867562, -0.885161821870766, -0.885178418039926, 
    -0.885187552637944, -0.885271540989251, -0.885266175732944, 
    -0.885686959032133, -0.885796426694203, -0.886732800166245, 
    41.7851467243973, 41.78520581764, 41.7854924854336, 41.7871110980829, 
    41.7871746095013, 41.7871983931199, 41.7871884798527, 41.7870973048765, 
    41.7870370230616, 41.7859610447487, 41.7859427714092, 41.785955093606, 
    41.7859354791523, 41.7859320425096, 41.7858909703717, 41.7858760497296, 
    41.7858794863084, 41.785599775486, 41.7854147818424, 41.7853933241996, 
    41.7852110974176, 41.7848115219888, 41.7851467243973), .Dim = c(23L, 
    2L))), class = c("XY", "POLYGON", "sfg"))), class = c("sfc_POLYGON", 
    "sfc"), precision = 0, bbox = structure(c(xmin = -0.886732800166245, 
    ymin = 41.4920576807899, xmax = -0.648597397154666, ymax = 41.7871983931199
    ), class = "bbox"), crs = structure(list(input = "WGS 84", 
        wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"latitude\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"longitude\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), row.names = 1:2, class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(itemID = NA_integer_, 
status = NA_integer_, Schema1 = NA_integer_, Schema2 = NA_integer_, 
Date = NA_integer_), .Label = c("constant", "aggregate", "identity"
), class = "factor"))

使用sf包很容易将此文件导出到 GML，只需：

st_write(mygeom, 'polygons.gml', driver='GML')

但是，当像这样导出时，我得到以下标准格式，它不符合我需要的标准。

<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
     xsi:schemaLocation="http://ogr.maptools.org/ polygons.xsd"
     xmlns:ogr="http://ogr.maptools.org/"
     xmlns:gml="http://www.opengis.net/gml">
  <gml:boundedBy>
    <gml:Box>
      <gml:coord><gml:X>-0.886732800166245</gml:X><gml:Y>41.4920576807899</gml:Y></gml:coord>
      <gml:coord><gml:X>-0.648597397154666</gml:X><gml:Y>41.7871983931199</gml:Y></gml:coord>
    </gml:Box>
  </gml:boundedBy>
                                                                                          
  <gml:featureMember>
    <ogr:polygons fid="polygons.0">
      <ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>-0.649062881756964,41.4921932026922 -0.649112159426678,41.4922493593042 -0.649217037157584,41.4923567934804 -0.649329802261076,41.4924515944872 -0.649385577124541,41.492509205028 -0.64940943118353,41.4925254105761 -0.649471971919662,41.4925678967614 -0.649418938777326,41.4926030387555 -0.649353483653088,41.4926842652395 -0.64927763496891,41.4927413468292 -0.64916144166702,41.4927247423209 -0.649060881454067,41.4926502665709 -0.648865049171286,41.4924908972756 -0.648768067415234,41.4924318736095 -0.648597397154666,41.4923280025819 -0.648913174239608,41.4920576807899 -0.648976607280396,41.4921104610734 -0.649027256609277,41.4921526044294 -0.649058695235838,41.4921884317747 -0.649062881756964,41.4921932026922</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
      <ogr:itemID>a1234</ogr:itemID>
      <ogr:status>1</ogr:status>
      <ogr:Schema1>A</ogr:Schema1>
      <ogr:Schema2>B</ogr:Schema2>
      <ogr:Date>2020-03-01</ogr:Date>
    </ogr:polygons>
  </gml:featureMember>
  <gml:featureMember>
    <ogr:polygons fid="polygons.1">
      <ogr:geometryProperty><gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>-0.886732800166245,41.7851467243973 -0.886458452915082,41.78520581764 -0.886077569154315,41.7854924854336 -0.885592339912652,41.7871110980829 -0.885526044922013,41.7871746095013 -0.885468424362897,41.7871983931199 -0.885411911891854,41.7871884798527 -0.885182818538855,41.7870973048765 -0.8850352039955,41.7870370230616 -0.885175736668516,41.7859610447487 -0.885147824091142,41.7859427714092 -0.885184202391753,41.785955093606 -0.885184872943689,41.7859354791523 -0.885184956762793,41.7859320425096 -0.885186297867562,41.7858909703717 -0.885161821870766,41.7858760497296 -0.885178418039926,41.7858794863084 -0.885187552637944,41.785599775486 -0.885271540989251,41.7854147818424 -0.885266175732944,41.7853933241996 -0.885686959032133,41.7852110974176 -0.885796426694203,41.7848115219888 -0.886732800166245,41.7851467243973</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></ogr:geometryProperty>
      <ogr:itemID>a27</ogr:itemID>
      <ogr:status>0</ogr:status>
      <ogr:Schema1>A</ogr:Schema1>
      <ogr:Schema2 xsi:nil="true"/>
      <ogr:Date>2020-03-01</ogr:Date>
    </ogr:polygons>
  </gml:featureMember>
</ogr:FeatureCollection>

在上面的文件中，上面的部分（1:12 行）是完全正确的，但其余的不是，因为我需要的是一个嵌套文件。我需要 Step 1 节点内的状态，以及 Step2 节点内的 Schemas。此外，当 Schema 2 为 NA 时，其节点消失。我可以使用 XML 包和下一个代码创建所需的格式。

library(XML) 
# create "GML"
mygml = newXMLDoc()
root = newXMLNode("featureCollection", parent = mygml)

for (i in 1:nrow(mygeom))
{
  
  # level 0
  l0<-newXMLNode("feature", attrs = c(fid=paste0('F',i)), parent = root)
  
  # level 1
  gp=newXMLNode('geometryProperty', parent = l0)
  s1 = newXMLNode('step1',  parent = l0)
  s2 = newXMLNode('step2',   parent = l0)
  
  # level 2
  gp1=newXMLNode('gml:Polygon',   attrs = c(srsName="EPSG:4326"),
                 namespaceDefinitions = c("gml" = "http://www.opengis.net/gml"),
                 parent = gp)
  
  s11=newXMLNode('item',  attrs=c(itemID=mygeom$itemID[i]), parent = s1)
  s21=newXMLNode('schema', parent = s2)
  s22=newXMLNode('schema', parent = s2)
  
  # level 3
  gp11=newXMLNode('gml:outerBoundaryIs', parent = gp1)
  
  s111=newXMLNode('Status',  mygeom$status[i], parent = s11)
  s112=newXMLNode('Date',  mygeom$Date[i], parent = s11)
  
  s211=newXMLNode('Schema',  mygeom$Schema1[i], parent = s21)
  s212=newXMLNode('Date',  mygeom$Date[i], parent = s21)
  
  if(!is.na(mygeom$Schema2[i])) {
    s221=newXMLNode('Schema',  mygeom$Schema2[i], parent = s22)
    s222=newXMLNode('Date',  mygeom$Date[i], parent = s22)
  }
  
  # level 4
  gp111=newXMLNode('gml:LinearRing', parent = gp11)
  
  # level 5 
  gp1111=newXMLNode('gml:coordinates', mygeom$geometry[i],  parent = gp111)
}

但是，这不是真正的 GML，即使在导出后我在这里编辑它也不会被识别。我的问题是，如何获得具有上述结构的 GML？

4

0 回答 0