Write Latitude and Longitude to Geotiff file

Question

I am basically trying to achieve the opposite of this question.

I have a set of latitude and longitude coordinates (with values) in the WGS84 coordinate system, that I would like to write to a geotiff (or just add to a gdal dataset) via the gdal python bindings.

For example, my starting data might be:

lat = np.array([45.345,56.267,23.425])
lon = np.array([134.689,128.774,111.956])
value = np.array([3.0,6.2,2.5])

How might one do this? Thanks!

are the locations of the data gridded or are they point locations on a map? — Mike T, Feb 11 '13 at 02:38

score 2 · Accepted Answer · edited Jul 24 '22 at 17:07

2

Although it is not in your question, it appears you need to project the lat/long data from the WGS84 datum to a UTM projection. This can be using the ogr2ogr command line from GDAL using the two options -a_srs 4326 -t_srs ???? (the target SRID). It can also be done internally with Python using the OGR module of GDAL. Here is an example of use.

There are two independent ways to get a raster from point data. The first is to interpolate the values in data, so that the values flood the region (or sometimes only the convex hull). There are many methods and tools to interpolate values in 2D. With GDAL, a comand-line tool gdal_grid is useful for this purpose, although I don't think it is possible to use from Python. Probably the simplest would be to use scipy.interpolate. Once you have a 2D NumPy array, it is simple to create a raster file with GDAL/Python.

The second method of converting the points to a raster is to burn the point locations to pixels on a raster. Unlike the first method, only the locations where the points are have values, while the values are not interpolated anywhere else in the raster. Rasterising or burning vectors into a raster can be done from a GDAL command line tool gdal_rasterize. It can also be done internally with GDAL/Python, here is an example.

edited Jul 24 '22 at 17:07

Glorfindel

21,988
13
81
109

answered Feb 11 '13 at 19:59

Mike T

41,085
18
152
203

I have found some code [HERE](http://monkut.webfactional.com/blog/archive/2012/5/2/understanding-raster-basic-gis-concepts-and-the-python-gdal-library/) that converts lat/lon points to UTM projected points and then pixels, which perhaps is doing some interpolation with `gdal.ApplyGeoTransform` to assign each point value across a given cell? – ryanjdillon Feb 12 '13 at 10:23
@shootingstars Do you just need a raster with a few pixels where the points exist, and empty space in-between? I.e., rasterize the points? What does the data represent? – Mike T Feb 12 '13 at 20:15
That's right. Once I've "rasterized" them, or have them in raster pixel units, I can use the inverse interpolation method I linked to in my first comment. Is this a bad way of going about things vs. just using the gdal_grid utility? – ryanjdillon Feb 13 '13 at 09:32
The data are whale sightings along ship transects cover the greater part of the North Atlantic. – ryanjdillon Feb 13 '13 at 13:55
@shootingstars Ok, I'm pretty sure you want to interpolate the data, not rasterise it (these are different). See updated answer. – Mike T Feb 13 '13 at 19:11

score 1 · Answer 2 · edited Jan 29 '15 at 17:36

It is possible to use gdal_grid from Python. I am using it. All you need to do is construct the command as if you were using it from the command line and put it inside a subprocess.call(com, shell=True). You need to import subprocess module first.

This is actually how I am using it:

pcall= "gdal_grid  --config 'NUM_THREADS=ALL_CPUS GDAL_CACHEMAX=2000'\
       -overwrite -a invdist:power=2.0:smoothing=2.0:radius1=360.0:radius2=360.0\
       -ot UInt16 -of GTiff -outsize %d %d -l %s -zfield 'Z' %s %s "%(npx, npy,\
       lname,ptshapefile,interprasterfile)

subprocess.call(pcall, shell= True)

NUM_THREADS option is available from gdal 1.10+

Write Latitude and Longitude to Geotiff file

2 Answers2