I have been working on this problem for a while now and have asked two seperate questions regarding it (Counting events from a GRIB raster, ValueError: array larger than output file, or offset off edge Python GDAL), that have helped me tremendously. Still, the final solution eludes me. The code should iterate through a 24 images and compare each value to the conditions set in the loop, then write the number of instances the conditions were met into a GeoTIFF image. Despite being as sure as I can be, that this code should work as intended, I have still not been able to produce anything more than a GeoTIFF full of zeroes.
import gdal
import numpy
import osr
import os
del_mapa = str(input('Vpiši pot do datotek v obliki "C:\\Mapa1\\Mapa2\\...\\Zadnja mapa": '))
os.chdir(del_mapa)
pix_vel = 1000
#Defines pixel size
ime_rez = str(input("Poimenuj datoteko z rezultati: ") + ".tif")
#User can provide a name for the output picture
format = "GTiff"
driver = gdal.GetDriverByName(format)
vir_mere = gdal.Open('inca_20140101-0100.grb')
trans = vir_mere.GetGeoTransform()
geotrans0 = trans[0]
geotrans1 = trans[1]
geotrans3 = trans[3]
geotrans5 = trans[5]
#Extracts transformation data from one of the GRIBs
dst_ds = driver.Create(ime_rez, 401, 301, 1, gdal.GDT_CInt32)
dst_ds.SetGeoTransform([geotrans0, geotrans1, 0, geotrans3, 0, geotrans5])
#and pastes it directly into the newly created GeoTIFF
projInfo = vir_mere.GetProjection()
srs = osr.SpatialReference()
srs.ImportFromWkt(vir_mere.GetProjectionRef())
dst_ds.SetProjection(srs.ExportToWkt())
raster = numpy.zeros((301, 401), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)
def krog1(baza,novi):
x = 0
while x < 25:
baza = gdal.Open(ime_rez)
if x < 10:
novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
#This should open the next file in line
if novi > 1 and novi < 12:
baza = baza + 1
else:
baza = baza + 0
elif x > 10:
novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
if novi > 1 and novi < 12:
baza = baza + 1
else:
baza = baza + 0
x += 1
return dst_ds.GetRasterBand(1).WriteArray(baza)
dst_ds = None
Edit: The second half with additions and corrections as per the advice provided.
raster = numpy.zeros((301, 401), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)
banda = vir_mere.GetRasterBand(1)
bandtip = gdal.GetDataTypeName(banda.DataType)
def krog1(baza,novi,bandtip):
vr_x = 0
while x < 25:
baza = gdal.Open(ime_rez)
baza2 = baza.GetRasterBand(1)
baza3 = baza2.ReadAsArray(0, 0, 301, 401).astype(numpy.CInt32)
if x < 10:
novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
#This should open the next file in line
re_raster = novi.GetRasterBand(1)
band_raster = re_raster.ReadAsArray(0, 0, 301, 401).astype(numpy.bandtip)
#Reads the raster as an array and then follows the condition checking
if band_raster > numpy.bandtip(1) and band_raster < numpy.bandtip(12):
baza = baza + 1
else:
baza = baza + 0
elif x > 10:
novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
re_raster = novi.GetRasterBand(1)
band_raster = re_raster.ReadAsArray(0, 0, 301, 401).astype(numpy.bandtip)
if band_raster > numpy.bandtip(1) and band_raster < numpy.bandtip(12):
baza = baza + 1
else:
baza = baza + 0
vr_x += 1
return dst_ds.GetRasterBand(1).WriteArray(baza)
dst_ds = None
Answer
As I've alluded to in my comments, you need to refine your function in terms of its arguments and outputs and then make sure you call it. That said, there are also some issues with the logic of your function and your raster access. I've rewritten your code as follows, with some comments for explanation:
import gdal
import numpy
import osr
import os
from osgeo import gdal_array # Used later to map data types
del_mapa = str(input('Vpiši pot do datotek v obliki "C:\\Mapa1\\Mapa2\\...\\Zadnja mapa": '))
os.chdir(del_mapa)
# User can provide a name for the output picture
ime_rez = str(input("Poimenuj datoteko z rezultati: ") + ".tif")
# Open your dataset to get the projection and transformation info
vir_mere = gdal.Open('inca_20140101-0100.grb')
# Get GetTransform - No need to break down into geotrans0 etc as the whole
# object can be passed to SetGeoTransform later
trans = vir_mere.GetGeoTransform()
projInfo = vir_mere.GetProjection()
# Get band GDAL datatype and map to Numpy type
v_band = vir_mere.GetRasterBand(1)
gdal_band_type = v_band.DataType
np_band_type = gdal_array.GDALTypeCodeToNumericTypeCode(gdal_band_type) # Maps the GDAL data type to a Numpy type
# Close band and dataset
v_band = None
vir_mere = None
# Define function - won't do anything until it's called
def krog1():
# Start by creating an array of zeros
arr = numpy.zeros((301, 401), dtype=np_band_type)
# Loop through the files in your current directory
for f in os.listdir(os.curdir):
# Open the file and read to array
ds = gdal.Open(f)
band = ds.GetRasterBand(1)
# Read entire raster into array
b_arr = band.ReadAsArray(0, 0, 401, 301).astype(np_band_type)
# Loop over each cell in the array and, if between 1 and 12, add 1 to the corresponding value in arr
for r in range(301):
for c in range(401):
if 1 < b_arr[r,c] < 12:
arr[r,c] += 1
# return the updated array
return arr
# Call the above function storing the output in out_arr
out_arr = krog1()
# Now create output and apply projection and transform info
dst_ds = driver.Create(ime_rez, 401, 301, 1, gdal.GDT_CInt32)
dst_ds.SetGeoTransform(trans)
srs = osr.SpatialReference()
srs.ImportFromWkt(projInfo)
dst_ds.SetProjection(srs.ExportToWkt())
# Write out_arr to band
dst_ds.GetRasterBand(1).WriteArray(out_arr)
# Flush to disk and close file
dst_ds.FlushCache()
dst_ds = None
This isn't necessarily the most efficient way but should hopefully highlight the issues you were having and give you more of an understanding of raster processing with GDAL. Again, I recommend this resource.
I've included gdal_array.GDALTypeCodeToNumericTypeCode(gdal_band_type)
here as this will give you the equivalent Numpy data type for your GDAL type which you can then use when creating or reading to arrays. Note that it requires the import of gdal_array
.
As you are dealing with user input for your file paths you should check that a valid path has been provided. I think this is why you had that NoneType error when calling 'baza.GetRasterBand(1)' - if GDAL can't open the file it silently returns a NoneType.
Update
The call to ReadAsArray
had the win_xsize
and win_ysize
arguments the wrong way round: b_arr = band.ReadAsArray(0, 0, 301, 401).astype(np_band_type)
should have been b_arr = band.ReadAsArray(0, 0, 401, 301).astype(np_band_type)
. Numpy shapes are always of the form (numrows,numcols) but the ReadAsArray
function needs the order of numcols,numrows. I've updated the above.
Alternatively, if you don't pass any arguments the whole array will be read by default i.e. b_arr = band.ReadAsArray().astype(np_band_type)
No comments:
Post a Comment