[GIS] Python – How to Get lat lon information from .asc file

esri-ascii-rastergeolocationlatitude longitudenumpypython

I have a .asc climate data file at hand and want to explore it to extract some information.

Based on searching, some general read methods have been tried as, import it as python pandas or numpy structure to see its contents.

Here the result:

if read in using pandas:

df = pd.read_table(path + "\max1951_1.asc")
df
                  NCOLS   1068
0                 NROWS    510
1    XLLCORNER  -141.0000
2    YLLCORNER   41.00000
3    CELLSIZE   0.83333330E-01
4      NODATA_VALUE  -999.00
5                      -999.00
6                      -999.00

Using numpy:

ascii_grid = np.loadtxt(path + "\max1951_1.asc", skiprows=6)
ascii_grid
array([[-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.],
       ...,
       [-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.]])

I still doubt what is the best solution or packet working with this kind of ESRI ASCII Raster format, not sure if the above print already shows all the file's content or not. Besides, my question is, if I want to extract any specific gridded cell's lat and lon information, How can I do that?

if asc file naturally contain those information so I can simply use certain existed packet to extract that?
Or the only method is use coordinate to calculate and convert base on the CELLSIZE and X,YLLCORNER.

Best Answer

You've tagged your question with ESRI-ASCII, why not use ESRI tools?

You can get properties of a raster or individual cell values. If you want to call these tools in python simply scroll to the bottom of these pages to see how to call them?

Related Solutions

[GIS] Get lat/lon points from image pixels

I find the simplest way that it can be done is from a georeferenced image by origin and cell size, if you can determine the Cell (row,column) on the image:

X = origin X + (Column * Cell Width) Y = origin Y + (Row * Cell height)

Usually.. the world file for a png (pgw) will give the top left cell and cell size in a 6 parameter transformation, be aware the Y cell size (height) is usually negative so it's a plus not a minus. Check out this site for a brief explanation of the format.

Basically it's: Cell Width 0 0 Cell height (usually negative) Upper left X Upper left Y

I'm sure that there are other ways to do this and inbuilt functions but the maths on this one are so easy! If the image isn't georeferenced then try to do so, it's not hard here is a link on how to do it in QGIS.

[GIS] Get information from geotiff with lat/lon

I modified your code slightly to print more than one value for a very small (clipped) version of a real raster:

from osgeo import gdal,ogr
from osgeo.gdalconst import *
import struct
import sys

#lat = float(sys.argv[2])
#lon = float(sys.argv[3])

def pt2fmt(pt):
    fmttypes = {
        GDT_Byte: 'B',
        GDT_Int16: 'h',
        GDT_UInt16: 'H',
        GDT_Int32: 'i',
        GDT_UInt32: 'I',
        GDT_Float32: 'f',
        GDT_Float64: 'f'
        }
    return fmttypes.get(pt, 'x')

my_path = "/home/zeito/pyqgis_data/test_long_lat.tif"

#ds = gdal.Open(sys.argv[1], GA_ReadOnly)
ds = gdal.Open(my_path, GA_ReadOnly)

if ds is None:
    print 'Failed open file'
#    sys.exit(1)

transf = ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount #1
band = ds.GetRasterBand(1)
bandtype = gdal.GetDataTypeName(band.DataType) #Int16
driver = ds.GetDriver().LongName #'GeoTIFF'
transfInv = gdal.InvGeoTransform(transf)

#(3,4), (5,2), (6,8)
p1 = [(16.4851469606, 44.5051091348), 
      (16.4672665386, 44.4863346918),
      (16.5215038186, 44.4782885019)]

for item in p1:
    px = (item[0]-transf[0])/transf[1]
    py = (item[1]-transf[3])/transf[5]
    #px, py = gdal.ApplyGeoTransform(transfInv, lon, lat)
    structval = band.ReadRaster(int(px), int(py), 1, 1,  buf_type = band.DataType )
    fmt = pt2fmt(band.DataType)
    intval = struct.unpack(fmt , structval)
    print "value: ", intval[0]

After running the code at the Python Console of QGIS, it could be observed that it works well (see next image). I corroborated that because I have my own version of Value Tool Plugin (print row and column indices where first pixel is 0,0) and another one that uses indices where first pixel is (1,1). Values printed at Python Console were as expected in all three cases.

Related Question