2

I tried to do this code Array to Raster image but there was the following error

ERROR 4: Unable to open myraster.tif to obtain file list.
ERROR 4: Unable to open EPSG support file gcs.csv.
Try setting the GDAL_DATA environment variable to point to the
directory containing EPSG csv files.

this is the code

import numpy as np
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import matplotlib.pylab as plt

array = np.array(( (0.1, 0.2, 0.3, 0.4),
                   (0.2, 0.3, 0.4, 0.5),
                   (0.3, 0.4, 0.5, 0.6),
                   (0.4, 0.5, 0.6, 0.7),
                   (0.5, 0.6, 0.7, 0.8) ))
# My image array      
lat = np.array(( (10.0, 10.0, 10.0, 10.0),
                 ( 9.5,  9.5,  9.5,  9.5),
                 ( 9.0,  9.0,  9.0,  9.0),
                 ( 8.5,  8.5,  8.5,  8.5),
                 ( 8.0,  8.0,  8.0,  8.0) ))
lon = np.array(( (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5),
                 (20.0, 20.5, 21.0, 21.5) ))
# For each pixel I know it's latitude and longitude.
# As you'll see below you only really need the coordinates of
# one corner, and the resolution of the file.

xmin,ymin,xmax,ymax = [lon.min(),lat.min(),lon.max(),lat.max()]
nrows,ncols = np.shape(array)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)   
# That's (top left x, w-e pixel resolution, rotation (0 if North is up), 
#         top left y, rotation (0 if North is up), n-s pixel resolution)
# I don't know why rotation is in twice???

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(4326)                     # This one specifies WGS84 lat long.
                                             # Anyone know how to specify the 
                                             # IAU2000:49900 Mars encoding?
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system 
                                                   # to the file
output_raster.GetRasterBand(1).WriteArray(array)   # Writes my array to the raster
  • What is the code that you ran to receive that output? You have a link but I think it is unfair on potential answerers to be required to follow it to try and synthesize your question. – PolyGeo Nov 21 '16 at 13:33
  • i use the code that is at the link if you click at "Array to Raster image" i will put the code now – Μανώλης Παναγιωτάκης Nov 21 '16 at 13:38
  • Did you cd (os.chdir(path)) into the directory that myraster.tif is in? – artwork21 Nov 21 '16 at 13:49
  • I use canopy to do the code so i just run it. do i have to create the file first? – Μανώλης Παναγιωτάκης Nov 21 '16 at 13:51
3

In my system, your code works well. It only needs this line:

output_raster = None

at the bottom.

After running it at the Python Console of QGIS, this is the result visualized at Map Canvas:

enter image description here

Editing Note:

I copy/paste your code in a text file (gdal_error.py) and, there were modified/added the lines into red rectangles (see next image). As you can see, it works perfectly in a QGIS 2.14-0-Essen (Windows system).

enter image description here

  • what line did you change? – Μανώλης Παναγιωτάκης Nov 21 '16 at 14:54
  • I added that line at the bottom of code. – xunilk Nov 21 '16 at 15:24
  • 1
    To add more info: That is a common GDAL Python Gotcha. An object needs to be dereferenced to be written to disk – Kersten Nov 21 '16 at 16:14
  • i ran it as a script and python console on QGIS and this is the error Traceback (most recent call last): File "<input>", line 1, in <module> File "c:/users/manwlas/appdata/local/temp/tmpuy4zqn.py", line 37, in <module> output_raster.SetGeoTransform(geotransform) # Specify its coordinates AttributeError: 'NoneType' object has no attribute 'SetGeoTransform' i have 2.14 Qgis and i work from windows. Do you think this is the problem? – Μανώλης Παναγιωτάκης Nov 22 '16 at 9:59
  • @Μανώλης Παναγιωτάκης Please, see my editing note. – xunilk Nov 22 '16 at 11:02

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.