1

I have a data file that gives me a porosity value and the lat/long that it applies to as a box.

For example

Porosity=14.58
90.0    0.0
90.0    360.0
89.0    360.0
89.0    0.0

I need this data in a format that ArcGIS can read and display as a map, preferably as a geotiff. What do I need to do to put this in a format that Arc can understand?

I have access to ArcGIS 10.4, Python, and MATLAB.

  • 1
    How proficient are you with arcpy? I would think that this would more readily convert to shapefile which can be rasterized. – Michael Stimson Aug 22 '17 at 5:13
1

You can use the Python OGR script provided below to convert your observation into polygon shapefile. It is just for one porosity observation (as an example).

After that, you can use some polygon to raster tool to convert the shapefile into GeoTIFF.

from osgeo import ogr
from osgeo import osr

def main():

    # read the input
    #observation = open("porosities.txt")
    #l0 = observation.read() # 'Porosity=14.58'
    #l1 = observation.read() # '90.0   0.0'
    #l2 = observation.read() # '90.0   360.0'
    #l3 = observation.read() # '89.0   360.0'
    #l4 = observation.read() # '89.0   0.0'

    # after parsing, spliting, converting, you would get..
    porosity = 14.58
    topLeft = (0.0, 90.0)
    topRight = (360.0, 90.0)
    bottomRight = (360.0, 89.0)
    bottomLeft = (0.0, 89.0)

    #------------------------------------------------------

    #create the layer
    outDriver = ogr.GetDriverByName("ESRI Shapefile")
    outDs = outDriver.CreateDataSource("Porosities.shp")
    outSrs = osr.SpatialReference()
    outSrs.SetWellKnownGeogCS("WGS84")
    outLayer = outDs.CreateLayer("Porosities", outSrs, ogr.wkbPolygon)
    #create the fields
    #outLayer.CreateField(ogr.FieldDefn("id", ogr.oftInteger))
    outLayer.CreateField(ogr.FieldDefn("porosity", ogr.OFTReal))

    outLayerDefn = outLayer.GetLayerDefn()
    #finally, create the feature
    f = ogr.Feature(outLayerDefn)
    f.SetField("porosity", porosity)

    polyRing = ogr.Geometry(ogr.wkbLinearRing)
    polyRing.AddPoint(*topLeft)
    polyRing.AddPoint(*topRight)
    polyRing.AddPoint(*bottomRight)
    polyRing.AddPoint(*bottomLeft)
    polyRing.AddPoint(*topLeft)

    geom = ogr.Geometry(ogr.wkbPolygon)
    geom.AddGeometry(polyRing)

    f.SetGeometry(geom)
    outLayer.CreateFeature(f)

    outFeature = None
    outDs = None

if __name__ == "__main__":
    main()

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.