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– PolyGeo ♦2016年11月21日 13:33:41 +00:00Commented 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Μανώλης Παναγιωτάκης– Μανώλης Παναγιωτάκης2016年11月21日 13:38:35 +00:00Commented Nov 21, 2016 at 13:38
-
Did you cd (os.chdir(path)) into the directory that myraster.tif is in?artwork21– artwork212016年11月21日 13:49:16 +00:00Commented 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?Μανώλης Παναγιωτάκης– Μανώλης Παναγιωτάκης2016年11月21日 13:51:03 +00:00Commented Nov 21, 2016 at 13:51
1 Answer 1
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:
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).
-
what line did you change?Μανώλης Παναγιωτάκης– Μανώλης Παναγιωτάκης2016年11月21日 14:54:33 +00:00Commented Nov 21, 2016 at 14:54
-
I added that line at the bottom of code.xunilk– xunilk2016年11月21日 15:24:41 +00:00Commented Nov 21, 2016 at 15:24
-
1To add more info: That is a common GDAL Python Gotcha. An object needs to be dereferenced to be written to diskKersten– Kersten2016年11月21日 16:14:11 +00:00Commented 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?Μανώλης Παναγιωτάκης– Μανώλης Παναγιωτάκης2016年11月22日 09:59:52 +00:00Commented Nov 22, 2016 at 9:59
-
@Μανώλης Παναγιωτάκης Please, see my editing note.xunilk– xunilk2016年11月22日 11:02:13 +00:00Commented Nov 22, 2016 at 11:02