GML file converted to ESRI shapefile with ogr2ogr missing spatial reference information

Question

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?

Answer 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: