4

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
4
  • 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, 2016 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, 2016 at 13:38
  • Did you cd (os.chdir(path)) into the directory that myraster.tif is in?
    – artwork21
    Nov 21, 2016 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, 2016 at 13:51

1 Answer 1

5

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

5
  • what line did you change? Nov 21, 2016 at 14:54
  • I added that line at the bottom of code.
    – xunilk
    Nov 21, 2016 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, 2016 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, 2016 at 9:59
  • @Μανώλης Παναγιωτάκης Please, see my editing note.
    – xunilk
    Nov 22, 2016 at 11:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.