lutraconsulting / serval

QGIS plugin for setting raster pixel values
38 stars 11 forks source link

GeoTiff raster #39

Closed epzt closed 2 years ago

epzt commented 2 years ago

Hi, I would like to use Serval plugin to apply modifications on raster layers based on GeoTiff format. When I want to do so, Serval complain it could not make any modifs on such format. I wold like to know the supported raster format by Serval to apply changes. Regards

epzt commented 2 years ago

For more information, here is the code I use to create my raster files:

gdal.SetConfigOption('GTIFF_DIRECT_IO','YES') gdal.SetConfigOption('GTIFF_VIRTUAL_MEM_IO','IF_ENOUGH_RAM') ds = gdal.GetDriverByName("GTiff").Create(_varname, xsize=nx, ysize=ny, bands=1, eType=gdal.GDT_Float64) ds.SetGeoTransform([extent.xMinimum(), pxsize, 0, extent.yMinimum(), 0, pysize]) ds.SetProjection(self.CRS.toWkt()) ds.GetRasterBand(1).WriteArray(_array)

Creating such a way, this raster is not recognized by Serval plugin. I did check in Serval code and didn't found any reason this raster coud not be change by Serval plugin. I did check "check_gdal_driver_create_option" and repeated all the operations of that function and this should return True. I then checked "check_layer" function and should also return True. I can't figure out where is my mistake. Regards

erpas commented 2 years ago

Hi Poizot,

How do you load the raster into QGIS? Is it stored on a local drive? Can you share an example?

epzt commented 2 years ago

Hi, here after the code I use now with some modifs since:

 ncdata = nc.Dataset(THENETCDFFILE)
 _array = ncdata["h"][:]
 _varname = "h"
 _x = ncdata["lon_rho"][:]
 _y = ncdata["lat_rho"][:]

 extent = QgsRectangle()
 extent.setXMinimum(np.nanmin(_x[0,:]))
 extent.setYMinimum(np.nanmin(_y[:,0]))
 extent.setXMaximum(np.nanmax(_x[0,:]))
 extent.setYMaximum(np.nanmax(_y[:,0]))

 pxsize = np.nanmean(np.diff(_x[0,:]))
 pysize = np.nanmean(np.diff(_y[:,0]))
 nx, ny = _array.shape

 filepath = os.path.join(os.path.dirname(THENETCDFFILE),f"{_varname}.tif")

 driver = gdal.GetDriverByName("GTiff")
 ds = driver.Create(filepath, xsize=ny, ysize=nx, bands=1, eType=gdal.GDT_Float32,
                                           options=["INTERLEAVE=BAND"])
 ds.SetGeoTransform([extent.xMinimum(), pxsize, 0, extent.yMinimum(), 0, pysize])
 srs = osr.SpatialReference()
 srs.ImportFromEPSG(4326)
 srs = srs.ExportToWkt()
 ds.SetProjection(srs)
 band = ds.GetRasterBand(1)
 band.WriteArray(_array)
 band.GetStatistics(0,1)
 band.SetNoDataValue(-9999)
 band.FlushCache()
 driver = None
 srs = None
 ds = None
 band = None
 del ds
 del band
 del driver
 del srs

Still have same problem with Serval plugin which do not recognized the created tif file as valid. After some tests, the problem is risen by Serval in the file "raster_handler.py", line193: "res = self.provider.setEditable(True)", the res variable is False or None. I tried to find information in the setEditable() function, but with no success. Still stuck on that actually.

To answer your question, I launch the raster through classical ways in QGIS, with the file explorer of QGIS and/or the raster manager. A example of created file can be download here: https://filesender.renater.fr/?s=download&token=c1654943-b485-4e3d-8b9b-6d57cf520428

Regards

epzt commented 2 years ago

Hi, just for information, if iI create the tif with the following options for the create function, Serval is allowed to change values.

  driver.Create(filepath, xsize=ny, ysize=nx, bands=1, eType=gdal.GDT_Float32,
                                      options=["INTERLEAVE=BAND","PROFILE=BASELINE"])

Of course, as a consequence there is no georeference of the tif created. Perhaps it will help to understand what going on. Regards

erpas commented 2 years ago

Make sure your pysize has a negative value and your script should produce a raster that can be turned editable in QGIS.

epzt commented 2 years ago

Thanks, works now. True pysize has to be negative. I can use now Serval to make my changes on rasters . Thanks for this plugin ;) Regards