Latitude/Longitude of each pixel using python and gdal

I needed lat/long of each pixel of a GeoTiff file. I searched through internet to find a solution and ,thanks to Stack Overflow, found a piece of code that I modified. As usual, python and gdal were used. Just follow the instructions below.

First get the number of rows and columns of you image by using gdalinfo (although it can be automated as well but lets just go with it)

gdalinfo raster.tif

Now edit the python code below and set rows and columns according to your image (don’t forget to change the image name).

from osgeo import gdal

# Open tif file
ds = gdal.Open('raster.tif')

# GDAL affine transform parameters, According to gdal documentation xoff/yoff are image left corner, a/e are pixel wight/height and b/d is rotation and is zero if image is north up. 
xoff, a, b, yoff, d, e = ds.GetGeoTransform()

def pixel2coord(x, y):
 """Returns global coordinates from pixel x, y coords"""
 xp = a * x + b * y + xoff
 yp = d * x + e * y + yoff
 return(xp, yp)

# get columns and rows of your image from gdalinfo
rows = 36+1
colms = 34+1

if __name__ == "__main__":
 for row in  range(0,rows):
  for col in  range(0,colms): 
   print pixel2coord(col,row)

Finally, execute this python script using following command (provided that you have named your file script.py)

python script.py

Here is the output of my raster.tif

(29.25, 6.0)
(29.5, 6.0)
(29.75, 6.0)
(30.0, 6.0)
(30.25, 6.0)
(30.5, 6.0)
(30.75, 6.0)
(31.0, 6.0)
(31.25, 6.0)
(31.5, 6.0)
.
.
.

In my case, I converted output to csv and opened it in qGIS along with raster.tif (shown below).

for_fb2

and a closer look

for_fb

 

Advertisements

Gdal – a swiss knife

572px-GDALLogoColor.svg

I worked with Global Inventory Modeling and Mapping Studies (GIMMS) dataset few weeks back, it is based on AVHRR and goes back to 1981.  For analysis, I had to resample GIMMS dataset from 8 km to 5 km (the version of GIMMS I was working with had 8 km resolution). I was working with approx. 736 images. Since I was performing an other task in R, thus, decided to combine resampling with it. I wrote a bash script to call my R script and passed source image and desired resolution image. The process was very simple as shown below (full script available here).

.......
# target image folder
src1folder="path-to-files"

#desired resolution image
src2folder="folder-for-desired-image"
src2img="desired-image-name.tif" 

for myfile in *.tif
do
  # target path and image name
  target="target-image-folder-path" 
  
  Rscript /media/external/muhammad/myscripts/scale.r $src1folder $myfile $src2folder $src2img $target
  
  # R temp folder is being filled by temp files, we need to clean it as well. 
  rm /tmp/R_raster_geoscilab/*.gri
  rm /tmp/R_raster_geoscilab/*.grd

done
.......

The R script was executed inside the loop for each image in folder “path-to-files”. R script (available here) read desired and source images and performed resampling using projectRaster() method. There were couple of other tasks performed by my R script but resampling was major.  When I tested this approach,  a single image took 12 minutes (5 images per hour). WOW! my script will run for 6 days to finish this task. Yes, scripts worked but 6 days were too much. I found out that it is projectReater() function which is taking almost 11 minutes. After a quick discussion with a friend of mine (R guru and big fan of gdal + Open Source) I tried gdal utilities. Quickly, I wrote a new script utilizing gdalwarp utility (command below, full script is available here).

gdalwarp -tr 0.050000000000000 -0.050000000000000 $outputfile $output_newfile

Notice the values for -tr switch which dictates the output image to be in 5km pixel size.  When I tested this image, WOW it took only 5 seconds for a single image. Still R was required for filtering but now I will not run my scripts for 6 days. It is important to note here that I used the default method for resampling in gdalwarp. Below is a GIMMS 5km resolution image (with nodata for further processing).

GIMMS-2011-dec-a

Moral of story: Try to use other tools rather than focusing on one.

Installing rgdal package for R 3

Day by day, my usage of R is increasing, it is simple, less code yet powerful. Recently I used R to get pixel values from a stack of images and generated graphs. Today i had to upgrade a server to R3 and ran into installation related issues

rgdal package installation was causing problems, error message was

gdal-config
--configure-args='--with-gdal-config=/usr/local/bin/gdal-config'
with appropriate values for your installation.

ERROR: configuration failed for package ‘rgdal’
* removing ‘/home/mycomputer/R/x86_64-pc-linux-gnu-library/3.0/rgdal’

The downloaded source packages are in
‘/tmp/Rtmp27M68g/downloaded_packages’

Warning message:
In install.packages("rgdal") :installation of package ‘rgdal’ had non-zero exit status

I didn’t have gdal-config on my computer and when i tried to install gdal-config following using command

sudo apt-get install libgdal1-dev

I got following errors

Reading package lists... Done
Building dependency tree 
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:

The following packages have unmet dependencies:
libgdal1-dev : Depends: libdap-dev but it is not going to be installed
E: Unable to correct problems, you have held broken packages.

I checked my gdal version using synaptic (alternate gdalinfo – -version command) and it was 1.9 (and this was the issue, gdal 1.10 was required) thus added ppa from ubutnugis.  Started synaptic package manager and clicked ‘Repositories’ from settings menu, went to ‘Other Software ‘ tab and added following repository (one by one)

deb http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu precise main

deb http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu precise main (Source Code)

Now from Edit menu, clicked ‘Reload Package Information’ , gdal 1.9 required upgrade, once update was complete  the command below worked

sudo apt-get install libgdal1-dev

relax it will take some time. Finally i installed package rgdal on R using

install.packages('rgdal')  

I checked installation by executing my graph R script and it worked.

rain-fall

But I discovered an other issue (after some time), now an applications (geonode) requiring gdal is not working. I installed geonode on an other server and this time instead of the above process, I did not update gdal but simply used the following command (use it after installing R)

sudo apt-get install libgdal1-dev libgdal-dev

Finally, I installed rgdal package. Till now I have not seen any issues with applications (and gdal version has not changed).