2

I am trying to convert Sentinel-2 cloud masks from GML format to ESRI shapefiles. I have managed to use ogr2ogr in a python script as described in the answer to Using ogr2ogr to convert GML to shapefile in Python? to do this. I downloaded ogr2ogr from http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ogr2ogr.py and my Python call to that looks like this:

ogr2ogr.main(["", "-f", "ESRI Shapefile", 'output_shapefile.shp', 'cloud_mask.gml'])

However when I load the new shapefiles into ArcMap, it tells me that they are "missing spatial reference information". I've opened the GMLs in a text editor and I can see the EPSG code.

GDAL isn't creating a .prj file at all. If I try with '-lco SPATIAL_INDEX=YES' it creates a .qix file (which it doesn't otherwise), but the problem remains. The output of ogrinfo -al -so for the GML file is as follows:

INFO: Open of 'cloud_mask.gml' using driver `GML' successful.
Metadata: NAME=MSK_CLOUDS pixels mask from data-strip S2A_OPER_MSK_CLOUDS_SGS__20190403T 205739_A019743_T13RGP_B00_MSIL1C

Layer name: MaskFeature
Geometry: Polygon
Feature Count: 4
Extent: (806760.000000, 3335340.000000) - (809760.000000, 3381900.000000)
Layer SRS WKT: (unknown)
gml_id: String (0.0) NOT NULL
maskType: String (6.0)

Is there some way I can get ogr2ogr to project the shapefiles that it creates correctly, or another way I can project them afterwards using some extra Python code?

3
  • did you try with -lco SPATIAL_INDEX=YES
    – nmtoken
    Jun 3 '19 at 16:02
  • What does ogrinfo -al -so report? How does the .prj file that GDAL creates look like?
    – user30184
    Jun 3 '19 at 16:51
  • @nmtoken and user30184 thanks for the comments - I've updated the question with the extra info you both asked for Jun 4 '19 at 7:54
1

I've had a similar issue, the only difference was that I wanted the Sentinel-3 footprint that I extracted from the xfdumanifest.xml. Fixed it by including a srsName tag with a CRS that was GML3+ compliant to my GML file, I'll paste it here for you if it can be of any help:

<gml:Polygon xmlns:gml="http://www.opengis.net/gml" srsName="http://www.opengis.net/def/crs/EPSG/0/4326">
  <gml:exterior>
    <gml:LinearRing>
      <gml:posList>0.303954 -69.9871 -0.209904 -70.1023 -0.723676 -70.2183 -1.23747 -70.3349 -1.75126 -70.4523 -2.26505 -70.5704 -2.77879 -70.6892 -3.2924 -70.8087 -3.80594 -70.9291 -4.31932 -71.0503 -4.83259 -71.1724 -5.34576 -71.2954 -5.85874 -71.4193 -6.3716 -71.5442 -6.69721 -71.6239 -7.13992 -69.8006 -7.5754 -67.9738 -8.00319 -66.1434 -8.42285 -64.3093 -8.83391 -62.4712 -9.23595 -60.629 -9.36789 -60.0141 -8.84874 -59.9015 -8.32955 -59.7888 -7.81021 -59.676 -7.29077 -59.5633 -6.77132 -59.4505 -6.2518 -59.3376 -5.73226 -59.2247 -5.21272 -59.1116 -4.69315 -58.9984 -4.17364 -58.885 -3.65426 -58.7714 -3.13501 -58.6575 -2.61583 -58.5434 -2.28612 -58.4709 -1.88064 -60.2908 -1.47325 -62.1099 -1.06438 -63.9282 -0.654437 -65.7461 -0.243832 -67.5637 0.16702 -69.3812 0.303954 -69.9871</gml:posList>
    </gml:LinearRing>
  </gml:exterior>
</gml:Polygon>

After applying your python call to ogr2ogr:

import ogr2ogr
ogr2ogr.main(["", "-f", "ESRI Shapefile", 'shp_test/footprint.shp', 'footprint.gml'])

I've ended up with this: enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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