Skip to main content
4 of 5
deleted 125 characters in body
Comrade Che
  • 8.2k
  • 32
  • 68

Also you could try to download txt files with footprint corners coordinates located here: https://ladsweb.modaps.eosdis.nasa.gov/archive/geoMeta/6/TERRA/

File example: https://ladsweb.modaps.eosdis.nasa.gov/archive/geoMeta/6/TERRA/2018/MOD03_2018-01-01.txt

*Python 3 Code to convert .txt MODIS footprints to ESRI shapefile (as @mikoontz requested). QGIS 3 required.

Disadvantage: 180 degree meridian problem does not solved.

enter image description here

import qgis.core
import qgis.PyQt.QtCore


def parse_txt_to_dict(txt_filepath):
    txt_as_dict = []
    with open(txt_filepath) as file:
        string = file.read()

    string = string.replace('\\n', '\n')
    lines_list = string.split('\n')

    # catch filed names
    fields_list = None
    for line in lines_list:
        if line.startswith("#") and 'GranuleID' in line:
            line = line.replace('# ', '')
            fields_list = line.split(',')
        # print(len(fields_list), fields_list)

    # catch values
    if fields_list is None:
        print('fields_list is None')
        return None

    for line in lines_list:
        if not line.startswith("#"):
            line = line.split(',')
            if not line == ['']:
                dictionary = dict(zip(fields_list, line))
                txt_as_dict.append(dictionary)

    return txt_as_dict


def dict_to_shp(dicts_list, shp_path):
    if dicts_list == []:
        print('Empty list. Skip.')
        return None

    # fields
    keys = dicts_list[0].keys()
    fields = qgis.core.QgsFields()
    for key in keys:
        fields.append(qgis.core.QgsField(key, qgis.PyQt.QtCore.QVariant.String))

    layer_name = 'test'
    # layer = qgis.core.QgsVectorLayer("Polygon?crs=epsg:4326", layer_name, "memory")
    qgis.core.QgsVectorLayer("Polygon?crs=epsg:4326", layer_name, "memory")

    writer = qgis.core.QgsVectorFileWriter(vectorFileName=shp_path,
                                           fileEncoding='UTF-8',
                                           fields=fields,
                                           geometryType=qgis.core.QgsWkbTypes().Polygon,
                                           srs=qgis.core.QgsCoordinateReferenceSystem("EPSG:4326"),
                                           driverName='ESRI Shapefile',
                                           )
    # geometry
    for dictionary in dicts_list:
        points = [qgis.core.QgsPointXY(float(dictionary['GRingLongitude1']), float(dictionary['GRingLatitude1'])),
                  qgis.core.QgsPointXY(float(dictionary['GRingLongitude2']), float(dictionary['GRingLatitude2'])),
                  qgis.core.QgsPointXY(float(dictionary['GRingLongitude3']), float(dictionary['GRingLatitude3'])),
                  qgis.core.QgsPointXY(float(dictionary['GRingLongitude4']), float(dictionary['GRingLatitude4'])),
                  qgis.core.QgsPointXY(float(dictionary['GRingLongitude1']), float(dictionary['GRingLatitude1'])),
                  ]

        feature = qgis.core.QgsFeature()
        feature.setGeometry(qgis.core.QgsGeometry.fromPolygonXY([points]))
        attr_vavalues = list(dictionary.values())
        feature.setAttributes(attr_vavalues)

        writer.addFeature(feature)
    del writer


def main():
    txt_path = r'C:\DELETE\GIS SE\txt\MOD03_2018-01-01.txt'
    shp_path = r'C:\DELETE\GIS SE\shp\MOD03_2018-01-01_shp.shp'
    dicts_list = parse_txt_to_dict(txt_path)
    dict_to_shp(dicts_list, shp_path=shp_path)


if __name__ == "__main__":
    main()
    print('Done!')
Comrade Che
  • 8.2k
  • 32
  • 68