0

I am working with a NetCDF file. One of the attributes of the file is its footprint given as just a list of GML coordinates e.g.

#ncds is a netCDF dataset
>>> ncds['METADATA/EOP_METADATA/om:featureOfInterest/eop:multiExtentOf/gml:surfaceMembers/gml:exterior']

<class 'netCDF4._netCDF4.Group'>
group /METADATA/EOP_METADATA/om:featureOfInterest/eop:multiExtentOf/gml:surfaceMembers/gml:exterior:
    gml:posList: 65.88228 88.01695 66.576965 91.344635 67.19473 94.85594 67.72871 98.5425 68.173035 102.3883 68.52232 106.3714 68.77156 110.4632 68.91724 114.6293 68.95755 118.83096 68.891426 123.02764 68.72009 127.17892 68.44618 131.24695 68.07342 135.19887 67.88231 136.86812 67.967384 136.9064 69.52042 137.66231 72.02265 139.18839 73.57038 140.42668 74.942116 141.82097 75.92515 143.06503 76.88609 144.55208 77.62853 145.93985 78.39927 147.66418 79.026024 149.33498 79.650375 151.30263 79.705925 151.49493 80.28056 153.68028 80.92471 156.64708 81.48355 159.81781 82.11964 164.40085 82.66919 169.65118 83.26404 177.84111 83.695564 -172.06956 83.90092 -155.65439 83.45412 -136.42804 82.56569 -122.955055 82.56569 -122.955055 83.802025 -116.99391 84.93152 -108.12669 85.865395 -94.649414 86.44607 -75.12895 86.4936 -51.877663 85.98876 -31.434368 85.09935 -17.029451 83.993744 -7.541537 82.77034 -1.1994907 81.47966 3.2010684 80.147964 6.36937 78.78965 8.722432 78.19511 9.556566 78.20421 9.9596405 78.267456 18.05271 77.957405 29.851692 77.536766 36.854866 77.03233 42.749462 76.60199 46.777275 76.12977 50.56255 75.73114 53.39411 75.28579 56.26016 74.89867 58.54604 74.48857 60.794594 74.4508 60.993713 74.04638 63.049828 73.558945 65.36215 73.09885 67.40047 72.51605 69.80719 71.9337 72.038734 71.143265 74.82091 70.2827 77.56101 68.971016 81.23353 67.285866 85.21143 65.88228 88.01695 65.88228 88.01695
    objectType: gml:LinearRing
    dimensions(sizes): 
    variables(dimensions): 
    groups:'

>>> footprint = [float(f) for f in getattr(ncds[r'METADATA/EOP_METADATA/om:featureOfInterest/eop:multiExtentOf/gml:surfaceMembers/gml:exterior'], 'gml:posList').split()]
>>> flon, flat = footprint[1::2], footprint[0::2] #get lats and lons 

I want to turn the list of coordinates into a valid WKT or PostGIS geometry like so:

MULTIPOLYGON (((88.01695 65.88228, 91.344635 66.576965, 94.85594 67.19473, 98.5425 67.72871, 102.3883 68.173035, 106.3714 68.52232, 110.4632 68.77156, 114.6293 68.91724, 118.83096 68.95755, 123.02764 68.891426, 127.17892 68.72009, 131.24695 68.44618, 135.19887 68.07342, 136.86812 67.88231, 136.9064 67.967384, 137.66231 69.52042, 139.18839 72.02265, 140.42668 73.57038, 141.82097 74.942116, 143.06503 75.92515, 144.55208 76.88609, 145.93985 77.62853, 147.66418 78.39927, 149.33498 79.026024, 151.30263 79.650375, 151.49493 79.705925, 153.68028 80.28056, 156.64708 80.92471, 159.81781 81.48355, 164.40085 82.11964, 169.65118 82.66919, 177.84111 83.26404, 180 83.35637644338722, 180 85.05115, 5.0188976408280670e-15 85.05115, -16.61581589418478 85.05115, -7.541537 83.993744, -1.1994907 82.77034, 4.2443893032536115e-15 82.4185304636445, 3.2010684 81.47966, 6.36937 80.147964, 8.722432 78.78965, 9.556566 78.19511, 9.9596405 78.20421, 18.05271 78.267456, 29.851692 77.957405, 36.854866 77.536766, 42.749462 77.03233, 46.777275 76.60199, 50.56255 76.12977, 53.39411 75.73114, 56.26016 75.28579, 58.54604 74.89867, 60.794594 74.48857, 60.993713 74.4508, 63.049828 74.04638, 65.36215 73.558945, 67.40047 73.09885, 69.80719 72.51605, 72.038734 71.9337, 74.82091 71.143265, 77.56101 70.2827, 81.23353 68.971016, 85.21143 67.285866, 88.01695 65.88228)), ((-122.955055 82.56569, -116.99391 83.802025, -108.12669 84.93152, -106.40024210506762 85.05115, -180 85.05115, -180 83.35637644338722, -172.06956 83.695564, -155.65439 83.90092, -136.42804 83.45412, -122.955055 82.56569)))

However my problem arises because they are not simple polygons. When plotted they self-intersect and actual fact are much better suited to being 2 separate polygons or a multipolygon.

How can I parse the data to give a valid polygon? For context this is in relation to the Sentinel-5P footprints, examples here. I was using this guide to try and help. The coordinate formatting in the first link is different to the netcdf file so the first code snippet is thing I am actually working with.

  • 1
    you need to convert the GML geometry into a geometry object and then serialize that object to WKT - rather than trying a search and replace op on the string of coordinates – Ian Turton Apr 30 at 7:38
  • i think my problem is that there is no valid GML geometry? In the file metadata it is literally just a string of coordinates – thesunnyscientist Apr 30 at 11:41
  • 1
    Then you are out of luck, file a bug report with whatever is generating the metadata – Ian Turton Apr 30 at 12:00

Your Answer

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

Browse other questions tagged or ask your own question.