[gdal-dev] GDAL slow to write GeoTIFF?

[jramm]

I'm writing a custom processing program using GDAL in C.

I'm processing a raster of roughly 150 000 * 200 000 pixels in windows of
256 * 256 pixels. 

I'm finding that traversing the raster and applying some basic processing to
each of the windows takes very little time; about 10 minutes on my machine.

However, as soon as start writing out each window (to the same output raster
file), the time goes up to about 30 hours

The writing is not special...just calling RasterIO like so:

GDALRasterIO(wt->hBand, GF_Write, xOff, yOff, xSize, ySize, data,
xSize, ySize, GDT_UInt32, 0, 0);


In a loop like so:

  while (movingwindow_next(&win) == 0)
{
process(winsize, &dg, &cats, &bt, resultBuffer);
xOff = movingwindow_xoff();
yOff = movingwindow_yoff();
//writer_writeWindow(&wt, resultBuffer, win.sizeX, win.sizeY, xOff,
yOff);
}

...I would've expected the code to slow down by 2, 3 or even 4 times, but
300 times is a bit crazy, no?

Are their any techniques I can use to optimise writing?

I am also not particularly limited by the output format if geotiff is not
the most efficient. I am only constrained by the fact that ESRI ArcGIS must
be able to read it...




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-slow-to-write-GeoTIFF-tp5203128.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL slow to write GeoTIFF?

[jramm]

By tiled do you mean split into N seperate TIFF files (in which case the
answer is no), or can a single TIFF contain a tiling scheme? If so, is this
some option I need to pass to the driver when creating the output?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-slow-to-write-GeoTIFF-tp5203128p5203202.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL slow to write GeoTIFF?

[jramm]

David, there should be no compression, but I will try setting it explicitly
in the driver...



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-slow-to-write-GeoTIFF-tp5203128p5203203.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Formats that support RAT's?

[jramm]

It seems GDAL will always let you create a Raster Attribute Table (in memory)
but may fail when writing to disk depending on whether the format supports
it? 

I can't seem to find any information anywhere on what formats support RAT's?
Is there a list I can refer to ??



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Formats-that-support-RAT-s-tp5203365.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Raster attribute tables in GDAL?

[jramm]

Does GDAL have any support for writing ESRI-style raster attribute tables?

These are the .vat.dbf files that often accompany geotiffs and other
formats




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Raster-attribute-tables-in-GDAL-tp5206653.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Raster attribute tables in GDAL?

[jramm]

That sounds great. Does this driver allow you to only write the dbf file
(and not the shp/shx etc files?)

I can't find much on the docs for doing that:
http://www.gdal.org/drv_shapefile.html

I'll have a look through the header files...

On 20 May 2015 at 15:38, Even Rouault-2 [via OSGeo.org] <
ml-node+s1560n5206664...@n6.nabble.com> wrote:

> Le mercredi 20 mai 2015 15:57:38, jramm a écrit :
> > Does GDAL have any support for writing ESRI-style raster attribute
> tables?
> >
> > These are the .vat.dbf files that often accompany geotiffs and other
> > formats
>
> I wasn't particularly aware of those sidecar files, but looking a bit I
> found
> one sample at https://dl.vecnet.org/files/j67313813
>
> It seems it is a standard .DBF file, so you can read/write it with the OGR
> Shapefile driver.
>
> Even
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> ___
> gdal-dev mailing list
> [hidden email] <http:///user/SendEmail.jtp?type=node&node=5206664&i=0>
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
> --
>  If you reply to this email, your message will be added to the discussion
> below:
>
> http://osgeo-org.1560.x6.nabble.com/Raster-attribute-tables-in-GDAL-tp5206653p5206664.html
>  To unsubscribe from Raster attribute tables in GDAL?, click here
> <http://osgeo-org.1560.x6.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=5206653&code=amFtZXNzcmFtbUBnbWFpbC5jb218NTIwNjY1M3wyMTIyMDM1OTQ5>
> .
> NAML
> <http://osgeo-org.1560.x6.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Raster-attribute-tables-in-GDAL-tp5206653p5206907.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL/OGR C API, how should I be dealing with 'handles'

[jramm]

Many C API functions return a typedef which is referred to as a handle. E.g.

OGRFeatureHOGRLayerHOGRGeometryH
These seem to occur when the C++ api returns a pointer.I have 3 questions:
Why doesnt the C API return pointers? 
How should I pass these 'handles' around. Should I create pointers to them?
Do I need to free them?




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-OGR-C-API-how-should-I-be-dealing-with-handles-tp5209726.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] gdal_rasterize does not burn polygon attributes

[jramm]

I'm running gdal_rasterize, trying to make it create a new tif from polygons,
using an attribute field as the burn value. Here is my command:

gdal_rasterize -tr 5 5 -a SOP -l test -a_nodata 0 test.shp test.tif
I end up with an image that is all 0's. The shapefile is correct and can be
read by OGR...what am I missing?




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-rasterize-does-not-burn-polygon-attributes-tp5210321.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL query creation options?

[jramm]

Is there a way to programatically query what dataset creation/layer creation
options are available for a driver?I'd like to be able to do something like:
driver->listDCO(*options)




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-query-creation-options-tp5210356.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] gdal_rasterize in multiple processes

[jramm]

I have a number of vector datasets to rasterize. 

I am finding that when the datasets are small, I have no problem starting
multiple gdal_rasterize processes (using python). 

If my datasets are large (I'm trying to create rasters with millions of
pixels), I'm finding that I can only start one gdal_rasterize process at
once, or sometimes, after a couple have started the rest refuse to start.

RAM shows usage below 50%, and CPU usage is also below 50%. 

I can't figure out whether this is a gdal_rasterize thing or not? 



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-rasterize-in-multiple-processes-tp5215591.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] gdal_rasterize increase memory usage

[jramm]

gdal_rasterize is limited to use just 10MB of memory (line ~640 of
gdalrasterize.cpp).

Is there anyway to change this (without having to recompile?)

I'm noticing that changing the output data format from Byte to Int16 or even
Int32 drastically reduces performance. This must be because the strict
memory limits imposed by gdal_rasterize means that the amount of data it
reads in one go is exponentially reduce each time I increase the integer
size. Thus it has to loop through the polygon set many more times. 

I'm seeing drastic slow down by just changing from Byte to Int16 when using
gdal_rasterize (8-10 times slower)

Given that even the cheapest desktop has at least 1GB of RAM, isn't the
scanline buffer size over conservative?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-rasterize-increase-memory-usage-tp5215862.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Where can I find ogr_FileGDB.dll v1.4 compilled for Windows?

[jramm]

You can get full binary distributions for windows (including) ogr_FileGDB.dll
from http://www.gisinternals.com/
They dont seem to have it in msvc_1800 releases, but in 1700 and 1600.

You will also need the FileGDB API dll from ESRI. You can get this from the
ESRI website (you will need to register, but it is free)



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Where-can-I-find-ogr-FileGDB-dll-v1-4-compilled-for-Windows-tp5220386p5220558.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] memory leak in gdal program

[jramm]

I'm using GDAL 201 to iterate features in a geodatabase feature class and add
a new text field, whose value depends on the integer value in another
field.It gives a fairly large memory leak (approx 1GB for every 1m features
iterated). I'm not sure where the leak isany ideas?The full code:int
main(int argc, char *argv[]){   // Get the input arguments  char *path =
argv[1];char *layer = argv[2];  char *field = argv[3];  // Open the 
dataset
in update mode  GDALAllRegister();  GDALDatasetH hDS;   hDS = 
GDALOpenEx(path,
GDAL_OF_VECTOR | GDAL_OF_UPDATE , NULL, NULL, NULL);if( hDS == NULL )   
{   
printf( "Open failed.\n" ); exit( 1 );  }   // Get the 
layerOGRLayerH hLayer;
hLayer = GDALDatasetGetLayerByName(hDS, layer); OGRFeatureDefnH hFDefn =
OGR_L_GetLayerDefn(hLayer); // Get the field index  int index =
OGR_FD_GetFieldIndex(hFDefn, field);// Create the new field OGRFieldDefnH
hFieldDefn; hFieldDefn = OGR_Fld_Create("depth_band", OFTString);
OGR_Fld_SetWidth(hFieldDefn, 32);if( OGR_L_CreateField( hLayer,
hFieldDefn, TRUE ) != OGRERR_NONE ){printf( "Creating depth_band
field failed.\n" );exit( 1 );}OGR_Fld_Destroy(hFieldDefn);
int value; // Field value   char depth_band[32]; //New attribute// Loop
through all featuresOGRFeatureH hFeature;   OGR_L_ResetReading(hLayer); 
int
count = 0;  while( (hFeature = OGR_L_GetNextFeature(hLayer)) != NULL)   
{   //
Get the depth band value and map it to the depth band description text  
value = OGR_F_GetFieldAsInteger(hFeature, index);   switch(value){  
case 1: 
depth_band = "x <= 0.3m";   break;  
case 2: depth_band = "0.3 < x <=
1m";break;  case 3: 
depth_band = "1 < x <= 3m"; break;  
case 4: 
depth_band = "3 < x <= 6m"; break;  
case 5: depth_band = "6 < x <=
9m";break;  case 6: 
depth_band = "x > 9m";  break;  }   
// Set the
new field value OGR_F_SetFieldString(hFeature,
OGR_F_GetFieldIndex(hFeature, "depth_band"), depth_band);   
count++;
printf("%d features processed\r", count);   }   GDALClose(hDS);
printf("\nProcessing complete!\n");}



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/memory-leak-in-gdal-program-tp5220559.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Transform raster based on values in shapefile attribute table

[jramm]

I have a raster (.tif), which has an accompanying shapefile attribute table
(.dbf). Basically, every value in the raster is a key into the attribute
table, which has a number of fields (there is a 'Value' field which is the
'key' corresponding to the raster values, then a number of attributes).

I want to derive a new raster based on a specific field in the attribute
table. I.e. for each pixel in the raster, find the corresponding row in the
attribute table, get the value of the desired field and replace the pixel
value with this new value. 

What (if anything) does GDAL provide to enable me to do this efficiently? Is
there an efficient method to look up rows (which in the OGR model I suppose
correspond to features) based on the value of a field?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Transform-raster-based-on-values-in-shapefile-attribute-table-tp5224659.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Proposed Block iterator API

[jramm]

Yes this would be incredibly useful. We typically process >100GB rasters and
so iterating in 'windows' is a must for us. 
We have typically implemented a solution outside of GDAL, with a 'next'
function that looks something like this:

int RasterWindow::next()
{
static int windowNoX = -1;
static int windowNoY = 0;

// Initialise the size of the window read
curSizeX = windowSize;
curSizeY = windowSize;

// Check which window number we are on
if (((windowNoY + 1) >= nY) && ((windowNoX + 1) >= nX))
{
// All windows have been read
return 1;
}
// We have reached the end of a row (x) of windows
else if ((windowNoX + 1) >= nX)
{
windowNoX = 0;
windowNoY += 1;
}
// There are still pixels left on the row...
else
{
windowNoX += 1;
}

// Calculate the pixel offsets
xOff = windowNoX * windowSize;
yOff = windowNoY * windowSize;

// Correct the current window size to make sure it won't run off the 
edge
of the raster
if (xOff + windowSize > RasterSizeX)
{
curSizeX = RasterSizeX - xOff;
}

if (yOff + windowSize > RasterSizeY)
{
curSizeY = RasterSizeY - yOff;
}

return 0;

}
It allows us to read more (or less) than the 'natural' block size (if
necessary) so we can send big lumps of data off for processing on a GPU.

For convenience and completeness, it would be helpful if the GetNext
function would optionally take a buffer which,  if present, would be
populated with the pixel values, negating the need for a line of boiler
plate to read the values. 

Another (possibly?) useful component could be some sort of pixel iterator.
For smaller/non-parallel processes, after reading a block we then loop over
the pixels, often with some boiler plate code to pull out the array indices
for various purposes:

uint32_t i, j, iGlobal, jGlobal;
for(pixelCount = 0 ; pixelCount < nPixels ; pixelCount++)
{
i = pixelCount % currentWindowSizeX;
j = (pixelCount - i) / currentWindowSizeX;
iGlobal = i + nXoff;
jGlobal = j + nYoff;

}

A further function which would take in the desired block size and handle
reading the blocks, allocating buffers etc behind the scenes and then simply
return individual pixels & indices would be quite interesting.

int GDALBlockIteratorGetNextPixel(GDALBlockIteratorH hIter, void *pixel, int
*i, int *j, int *iGlobal, int *jGlobal);



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-dev-Proposed-Block-iterator-API-tp5237743p5237882.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL Postgis Driver using SPI?

[jramm]

Hi

I am working on a project to manipulate huge rasters in a postgis database.
In order to achieve best performance, we are implementing most of our code
'server-side', e.g. as a postgresql extension. 

It occurs to me that another GDAL postgis driver using SPI
(http://www.postgresql.org/docs/current/static/spi.html) instead of libpq
would be entirely feasible and offer better performance as it executes
directly against the database. Making it perfect in situations where:

- The postgres database is local (so using libpq to setup a network
connection is not necessarily the best method)
- You want best read/write performance and are happy to only operate gdal
'server side'
- You want to use gdal in order to write postgres/postgis extensions in C
and expose them to clients as SQL (e.g. thin clients, webapps etc..)

Using SPI would require building the driver against postgres and is somewhat
similar to using libpq, such that much of the existing code in the postgis
raster driver can be reused. 

Any thoughts on this?




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-Postgis-Driver-using-SPI-tp5247264.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Read data without reading 'null data'

[jramm]

We are storing large raster datasets (~500,000 x 500,000) in geotiff files.
The array data is typically very sparse (density is 0.02% or less) and
compression greatly reduces the geotiff size.

However, when processing data we read the data in chunks, which are
automatically decompressed..this means we are often reading in large areas
of no data and 'skipping' them. 
Is there a way to only read data with values so we may apply the processing
to these cells, and write them back out, skipping the no data? I guess this
would be equivalent to just not decompressing?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Read-data-without-reading-null-data-tp5251626.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Standard GDAL format for HDF5?

[jramm]

I was wondering if there would be any interest in a standard GDAL format for
HDF5 files?

The major problem I have with current formats supported by the HDF5 driver
is that they are very application-specific, with requirements on the
structure of groups and group names. 

What I propose would make no assumptions about the file grouping structure,
but require that some attributes are set on the dataset itself:

- Geotransform
- Projection
- NoDataValue

To read a certain dataset, you would pass the driver the 'posix style' path
to the dataset (as per HDF5 docs). I'm not sure how this would fit with the
driver, but perhaps something similar that is done for postgres  (i.e. all
in one string then the elements parsed out):

"/My/Path/To/file.hdf5/path/to/dataset"

This would allow any file structure and for multiple datasets within a file
to support different projections, geotransforms etc etc

Any thoughts?




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Standard-GDAL-format-for-HDF5-tp5253334.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Implementing write support for GDAL driver

[jramm]

Hi
I am implementing a write driver for the Vertical Mapper/Northwood Grid
format (NWT_GRD), which already has read support. 
I have attached a diff of what I have got so far. Needless to say, it is not
complete or working :D.

I'm looking for information on development, specifically:

- Precisely what functions need to be implemented so the writing is fully
supported? (i.e. the minimum to be able to create a dataset from scratch and
write data to it using RasterIO). 
- Is there any info out there on debugging GDAL drivers? What is the
standard approach. Currently im trying to use 'gdal_translate' as a way in
to get it to do something so I can set some break points...perhaps writing a
little test program would be a better method. Which brings me to
- What standards are there for writing unit tests, test programs? I can see
a load of python under 'autotest' which im not familiar with. Are there any
guidelines on thise (the README is none too helpful).

Ive been over the driver implementation tutorial a number of times but info
on implementing write support is fairly scant. I'd happily add to this if I
can get this done..

Thanks for any help
grddataset_diff.txt
  



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Implementing-write-support-for-GDAL-driver-tp5258855.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Implementing write support for GDAL driver

[jramm]

Thanks a lot
Thats useful and encouraging!
I have not implemented CreateCopy -. Ive not looked to closely at
gdal_translate, but if CreateCopy is not available will it fall back on
Create and the write methods or is this a likely cause of errors? From
first impressions it looks like CreateCopy is not entirely necessary -
mostly a shortcut when writing a whole dataset to another, but I've not
looked closely at it yet.




On 30 March 2016 at 16:36, Even Rouault-2 [via OSGeo.org] <
ml-node+s1560n5258858...@n6.nabble.com> wrote:

> Le mercredi 30 mars 2016 17:17:11, jramm a écrit :
>
> > Hi
> > I am implementing a write driver for the Vertical Mapper/Northwood Grid
> > format (NWT_GRD), which already has read support.
> > I have attached a diff of what I have got so far. Needless to say, it is
> > not complete or working :D.
> >
> > I'm looking for information on development, specifically:
> >
> > - Precisely what functions need to be implemented so the writing is
> fully
> > supported? (i.e. the minimum to be able to create a dataset from scratch
> > and write data to it using RasterIO).
>
> What you implemented in your patch looks good : IWriteBlock,
> SetGeoTransform,
> SetProjectionRef, Create.
> With that, you can have full "dynamic creation" (as named in
> http://gdal.org/gdal_drivertut.html#gdal_drivertut_creation ), thus
> allowing
> it to be used as the target of gdalwarp for example.
>
> > - Is there any info out there on debugging GDAL drivers? What is the
> > standard approach. Currently im trying to use 'gdal_translate' as a way
> in
> > to get it to do something so I can set some break points...perhaps
> writing
> > a little test program would be a better method.
>
> Using gdal_translate is a reasonable way to test & debug write support.
> When
> used without any particular arguments, it's mostly a wrapper of
> CreateCopy().
>
> > Which brings me to
> > - What standards are there for writing unit tests, test programs? I can
> see
> > a load of python under 'autotest' which im not familiar with. Are there
> any
> > guidelines on thise (the README is none too helpful).
>
> Not really. Generally looking at existing tests should give you a good
> idea of
> how to extend. The basic is to add a method that returns 'success' or
> 'fail'
> strings, and to add it to the gdaltest_list array.
>
> The gdaltest class has GDALTest() and testCreateCopy() methods that make
> it
> generally easy to test write support. Have a look at
> autotest/gdrivers/bt.py
> for example.
>
> I guess you found the existing tests for NWG_GRD are in
> autotest/gdrivers/nwt_grd.py
>
> >
> > Ive been over the driver implementation tutorial a number of times but
> info
> > on implementing write support is fairly scant. I'd happily add to this
> if I
> > can get this done..
> >
> > Thanks for any help
> > grddataset_diff.txt
> > <http://osgeo-org.1560.x6.nabble.com/file/n5258855/grddataset_diff.txt>
> >
> >
> >
> > --
> > View this message in context:
> >
> http://osgeo-org.1560.x6.nabble.com/Implementing-write-support-for-GDAL-dr
> > iver-tp5258855.html Sent from the GDAL - Dev mailing list archive at
> > Nabble.com.
> > ___
> > gdal-dev mailing list
> > [hidden email] <http:///user/SendEmail.jtp?type=node&node=5258858&i=0>
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> ___
> gdal-dev mailing list
> [hidden email] <http:///user/SendEmail.jtp?type=node&node=5258858&i=1>
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> --
> If you reply to this email, your message will be added to the discussion
> below:
>
> http://osgeo-org.1560.x6.nabble.com/Implementing-write-support-for-GDAL-driver-tp5258855p5258858.html
> To unsubscribe from Implementing write support for GDAL driver, click here
> <http://osgeo-org.1560.x6.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=5258855&code=amFtZXNzcmFtbUBnbWFpbC5jb218NTI1ODg1NXwyMTIyMDM1OTQ5>
> .
> NAML
> <http://osgeo-org.1560.x6.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Implementing-write-support-for-GDAL-driver-tp5258855p5258885.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] How to consider data types when implementing IWriteBlock

[jramm]

Hi
I am adding write support for the Northwood Grid (Vertical Mapper/MapInfo
format) driver. 
This format only allows 32 bit float as the datatype, and to make things
more tricky it actually stores this on disk as either 16 or 32 bit ints,
using scaling rules.

The current read driver applies the scaling to each pixel when reading, to
return a buffer of 32 bit floats and I look to do the reverse when writing. 

When implementing `WriteBlock`, is it OK to assume that pImage will always
be of the correct type? (float 32). I've noticed in debugging using
gdal_translate that if I dont specify the output data type, it carries on
and passes the source buffer (which is not float 32) to IWriteBlock of the
destination. 

Given that only allowing 32bit floats is quite restrictive, should I:

- raise an error in the `Create` method if the data type isnt GDT_Float32
- Attempt to handle all gdal data type in IWriteBlock by storing the user
requested data type in some variable (poDS->eUserDataType) and then using
GDALCopyWords within IWriteBlock to first copy pImage to a buffer of 32 bit
floats before proceeding?

Thanks




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/How-to-consider-data-types-when-implementing-IWriteBlock-tp5260901.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] How to consider data types when implementing IWriteBlock

[jramm]

Hi
I have implemented the write support in the NWT_GRD format, having to make a
few compromises. I'd like to go through them and see if they are reasonable.
I have also attached the diff.

1. The existing read driver reports 4 bands for a GRD dataset. The first 3
are essentially 'virtual' - they are RGB bands which are interpreted from
the single band of data actually on disk using color info in the header.
When opening in write mode, I have prevented more than 1 band being created.
Therefore in update mode, you cannot read from the 'virtual' bands; it is
created with 1 band to which you write the actual data. So you have some
slightly odd behaviour in that, you open it in `update` mode, it has 1 band,
you then close and reopen in `read only` mode, it has 4 bands. I wonder if
it would be better for the read driver to utilize the GDAL ColorTable stuff
somehow, instead of creating these 'virtual' bands?

2. The no data value is fixed on disk as 0; this is then converted to
-1.e37f by the read driver (as per the vertical mapper grid standard) upon
reading. On writing, I have allowed setting the no data value to anything,
and these values are treated as no data when writing. So you may create a
dataset having set the no data value to e.g. -999. However, when you close
it and open it in 'read' mode, the nodata value will be reported as -1.e37f

3. The GRD format supports data being written on disk as either 16bit
unsigned ints or 32bit unsigned ints, In the NWT_GRDRasterBand constructor,
allowance is made for this, however in the read driver, values are always
read as unsigned short (suggesting that if data were 32bit ints, we might
encounter over/underflow). I have done the same in the IWriteBlock, but I am
wondering if their should be a check to see whether the written value should
be an `unsigned short` or an `unsigned int`

4. Spatial Reference. The spatial reference is stored on disk in the 'MITAB'
format and a conversion function is used in the existing GetProjectionRef. A
function to convert the other way exists and I use this in SetProjection.
However, some stuff is obviously lost in the conversion as setting the
projection with e.g. EPSG:27700 and then reading it back loses the EPSG
information. Is there a way (using PamDataset maybe?) to store the GDAL
projection so that GDAL will always recognise it correctly, while still
writing the NWT_GRD format projection in the correct place in the tab/grd
file?

Are these reasonable compromises or can it be done better?

Thanks
James



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/How-to-consider-data-types-when-implementing-IWriteBlock-tp5260901p5261675.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] How to consider data types when implementing IWriteBlock

[jramm]

Even Rouault-2 wrote
> Skimming through the code, it seems you must provide zmin / zmax to do the 
> scaling of floating point elevations to integral values. So either it is
> user 
> provided, or you do a prior statistics computation in CreateCopy() if the 
> input dataset hasn't min/max metadata.
> In any case, checking overflow and warning if it occurs might be a good
> idea.

So far, I haven't implemented CreateCopy, although this would probably be a
very good idea to allow min/max to be set automatically & properly.
In Create(), it just defaults the min/max to SHRT_MIN and SHRT_MAX. Would
the creation options be a good way to expose it to the user, to be set
explicitly?


Even Rouault-2 wrote
> Yes using PAM might be an option. Basically our SetProjection() should
> also 
> call GDALPamDataset::SetProjection(). And on reading you first check if 
> GDALPamDataset::GetProjectionRef() returns a non-empty string, in which
> case 
> you return it. Otherwise you return the MITAB decoded one. 

Great. I will try doing this. 



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/How-to-consider-data-types-when-implementing-IWriteBlock-tp5260901p5261827.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Geotiff FillEmptyTiles with no data value?

[jramm]

Hi

When writing geotiffs, if I dont write blocks they will automatically be
filled on close by the FillEmptyTiles.
It appears that this will only fill with zeros - is it possible to make it
fill with the no data value instead?

This is potentially a huge time saver when processing a large, fairly sparse
geotiff using python as relying on gdal to fill in 'null' blocks appears far
quicker than calling the python write method. 




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Geotiff-FillEmptyTiles-with-no-data-value-tp5264310.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Segfault in GDALWriteBlock

[jramm]

Hi I am getting a segfault in a call to GDALWriteBlock with this code:



Stepping through seems to suggest that the segfault is arising from line 171
of "geo_new.c" in libgeotiff:



I cant find previsely what it is I have done wrong to cause this in my
codeany ideas?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Segfault-in-GDALWriteBlock-tp5264660.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Segfault in GDALWriteBlock

[jramm]

im building off trunk, so Gdal2.1dev

Ill get that stack trace...



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Segfault-in-GDALWriteBlock-tp5264660p5264899.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Geotiff FillEmptyTiles with no data value?

[jramm]

NoDataValue is set on a band level and FillEmptyTiles operates at a dataset
level. 

I've never heard of a geotiff with different nodata per band - is this even
possible, or would it be reasonable to take the nodatavalue of the 1st band
and use this to initialise the empty tiles block?

Then there is the problem of the block array being initialised as GBYte:

GByte *pabyData = (GByte *) VSI_CALLOC_VERBOSE(nBlockBytes,1);

NoDataValue is always returned as a double, but actual data type may be
different values...Am I right in thinking TIFFTileSize will return the
number of bytes in a block for the data type of the tiff file (e.g. so block
size = nBlockBytes / sizeof(data_type)?

Would something like this be a reasonable way to initialise pabyData

int typeSize;
switch (eType){
  case GDT_Byte:
  typeSize = sizeof(GByte)
  case GDT_Int16:
  typeSize = sizeof(short)

   ...
}


for (int i =0; i < nBlockBytes; i+=typeSize)
{
*pabyData = nullValue;
}



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Geotiff-FillEmptyTiles-with-no-data-value-tp5264310p5264905.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Geotiff FillEmptyTiles with no data value?

[jramm]

Hi
I have made changes to FillEmptyTiles so that if nodata is set, then it will
always fill with nodata, otherwise 0. 

I have attached the raw diff...I have no idea how to submit a change
request/review etc?

fillempty_nodata.diff
  



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Geotiff-FillEmptyTiles-with-no-data-value-tp5264310p5265514.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] How to consider data types when implementing IWriteBlock

[jramm]

As an update to this; I have finished write support for the GRD driver for
now, with ZMAX/ZMIN exposed as creation options and will be explicitly
calculated in CreateCopy if not passed in - exactly as you suggested Even,
thanks. 

I'm waiting permission to get a osgeo user id so I can raise a ticket and
submit the patch :D



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/How-to-consider-data-types-when-implementing-IWriteBlock-tp5260901p5268160.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Read data without reading 'null data'

[jramm]

Hi 
I am having a few problems with using this (via python).

I have a test dataset with a 16 x 16 blocksize. In the below examples, bnd =
ds.GetRasterBand(1)

First, I try getting the blocksizes using the metadata:



This gives:



None of the first 10 blocks are empty - this is strange as it is a
handwritten dataset and i made sure there were some empty blocks.

So, I try reading each block and then getting the number of elements that
are not equal to the no data value, using numpy:


This gives:


blocks 9 and 10 are empty, furthermore, the non empty sizes are completely
different...what size is the metadata reporting here?
Am I calculating offsets wrong?
Even if that is the case, a full block should never be larger than 256. 
returns 




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Read-data-without-reading-null-data-tp5251626p5270080.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Read data without reading 'null data'

[jramm]

This is my mistake. 
I had not truly created a sparse raster! I keep forgetting that the gdal
tools (gdal_translate) will still write the sparse blocks and the sparse
creation option is just to indicate that FillEmptyBlocks should be skipped. 
Creating the sparse file by hand correctly sets the block metadata. 

I can see that BLOCK_SIZE_{0}_{1} will return the size in bytes, or null.
What is interesting then is that if the raster isnt sparse (As in the first
example), BLOCK_SIZE_{0}_{1} appears to be returning garbage...



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Read-data-without-reading-null-data-tp5251626p5270084.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Discover whether a GDALDataset is raster or vector?

[jramm]

Hi
Given a GDALDataset pointer, which has been opened using GDALOpenEx(...),
what is the best way to discover whether the dataset is raster or vector?

I have thought of checking the drivers' metadata for GDAL_DCAP_OF_RASTER or
VECTOR, but this could potentially return YES for both if it is e.g the
geopackage driver?





--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Discover-whether-a-GDALDataset-is-raster-or-vector-tp5270223.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Pass clipsrc to VectorTranslateOptions

[jramm]

Im using VectorTranslate and VectorTranslateOptions in python bindings.

It seems like there is nothing equivalent to the ogr2ogr -clipsrc option in
VectorTranslateOptions.

Is there another way of passing this in?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Pass-clipsrc-to-VectorTranslateOptions-tp5271367.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Create standalone PAM datasets?

[jramm]

Does GDAL offer any functionality (in the python bindings) to create a PAM
dataset (.aux.xml file) without a 'main' dataset attached?

This is a useful metadata container that we would like to use for other
formats (e.g. non spatial CSV file, which might share the same custom
metadata as a geotiff) we derive from spatial data.





--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Create-standalone-PAM-datasets-tp5272004.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdalbuildvrt for multiband files?

[jramm]

Two options:

1. Use the xml module in python to build up the VRT file. You would create a
VRTRasterBand node for each band of each file. Note you will need to figure
out the maximum extent and raster size (in pixels) first, so this will be
easiest if you enforce the same pixel size and projection on your inputs

2. Create using the GDAL library. Something along the lines of

   drv = gdal.GetDriverByName("VRT")
   dst =  drv.Create()
   for src in src_list:
   for i in range(src.RasterCount):
  dst.AddBand(...)

Again, you would need to work out the size of the final raster beforehand.

The start of a very naive python implementation using xml is below. Look at
http://www.gdal.org/gdal_vrttut.html
to see what further nodes and attributes need creating 

from xml.etree.ElementTree import *
def stack(src_list):
# Get the maximum extent
ulx = []
uly = []
lrx = []
lry = []
xres = []
yres = []
for src in src_list:
   gt = src.GetGeoTransform()
   ulx.append(gt[0])
   uly.append(gt[3])
   lrx.append(gt[0] + (gt[1] * src.RasterXSize))
   lry.append(gt[3] + (gt[5] * src.RasterYSize))
   xres.append(xres)
   yres.append(yres)
   # Dont handle skewed images
   assert gt[2] == gt[4] == 0

   # Check resolutions are all the same
   assert xres.count(xres[0]) == len(xres)
   assert yres.count(yres[0]) == len(yres)

   dst_ulx = min(ulx)
   dst_lrx = max(lrx)
   dst_uly = max(uly) # for north up image
   dst_lry = min(lry)
   nxpixels = int((dst_lrx - dst_ulx) / xres[0])
   nypixels = int((dst_lry - dst_uly) / yres[0])

root = Element("VRTDataset")
root.set("RasterXSize", str(nxpixels))
root.set("RasterYSize", str(nypixels))

gt_node = SubElement(root, "GeoTransform")
gt_node.text = ", ".join([dst_ulx, xres[0], 0, dst_uly, 0, yres[0]])

for src in src_list:
   for i in range(src.RasterCount):
   bnd_node = SubElement(root, "VRTRasterBand")
   src_node = SubElement(bnd_node, "SimpleSource")
   pth_node = SubElement(src_node, "SourceFilename")
   pth_node.text = src.GetDescription()
   srcbnd_node = SubElement(src_node, "SourceBand")
   srcbnd_node.text = str(i+1)
  
A possibly simpler option would be to use gdalbuildvrt to create a VRT for
the first band of each dataset, then read the VRT using python and duplicate
the relevant nodes to add in new bands, then write it out again.





--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-dev-gdalbuildvrt-for-multiband-files-tp5272094p5272472.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] File Geodatabase raster driver

[jramm]

Based on Even's work on FGDB vector, I have begun looking at raster data,
here: https://github.com/JamesRamm/fgdb_raster_spec

So far everything I've found is pretty much based on the output of v the
dump_gdb script. I have two troubling areas currently blocking progress 

1. I can locate pixel data blocks but cannot decipher how to read them, even
when the contents are known.

2. It is unclear how more than one raster is handled. It appears to be based
solely on the names of files but any patterns/groupings are not obvious 

Does anyone have any more experience with this or any tips on how to
proceed?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/File-Geodatabase-raster-driver-tp5275154.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Call for discussion on RFC 63 : Sparse datasets improvements

[jramm]

Even Rouault-2 wrote
> For GeoTIFF, the unit of query remains the block. So if a block is missing
> is 
> present, it computes how many pixels of the window of request intersect
> the 
> block, and count them as present.
> 
> Imagine that you have a raster of dimensions 20x20, with tiles of size
> 10x10. 
> Imagine the first block (pixels 0,0 to 9,9) is missing. And you call 
> GetDataCoverageStatus() with (5,5,15,15). Then you will have a
> contribution of 
> 10x5 valid pixels from the top-right block, 5x10 from the bottom-left
> block 
> and 10x10 from the bottom-right block, hence a non-null percentage of 200.
> / 
> 225. * 100 = 88.9 %

Ah ok, makes sense. 
Having SPARSE_OK=YES mean that gdal_translate will actually skip null blocks
is a good additiona also - more intuitive. 





--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-dev-Call-for-discussion-on-RFC-63-Sparse-datasets-improvements-tp5275443p5275600.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] GDAL raster processing library?

[jramm]

I wonder if there is a use case/interest in a small library within GDAL to
facilitate generic raster processing. 

My idea would be to have a user-extensible framework to run processing
functions on whole rasters, with a growing library of common-such functions
within GDAL. 

Something along the lines of this:

/***/
typedef CPLErr (*GDALRasterProcessFn)(double *padfInArray, double
*padfOutArray, 
int nWindowXSize, int nWindowYSize, void *pData, double
*pdfNoDataValue);
/**
* \brief Definition of a raster processing function. 
*
* A GDALRasterProcessFn accepts an array of data as input, applies custom
logic and writes the output to padfOutArray.
* Such a function can be passed to GDALRunRasterProcess to apply custom
processing to a GDALDataset in chunks and create 
* a new GDALDataset. 
*
* @param padfInArray The input array of data. 
*
* @param padfOutArray The output array of data. On first call (via
GDALRunRasterProcess) this will be an empty, initialised array,
*which should be populated with the result of calculations on
padfInArray. In subsequent calls it will contain the result of the 
*previous window. 
*
* @param nWindowXSize the actual x size (width) of the read window.
*
* @param nWindowYSize the actual y size (height) of the read window. The
length of padfInArray == padfOutArray == nWindowXSize * nWindowYSize
*
* @param pData Process-specific data. This data is passed straight through
to the GDALRasterProcessFn and may contain e.g user defined parameters.
* The GDALRasterProcessFn definition would define the structure/type of
such data.
*
* @param pdfNoDataValue The no data value of the dataset
*/

CPLErr GDALRunRasterProcess(GDALRasterProcessFn processFn, GDALDataset
*poSrcDataset, 
GDALDataset *poDstDataset, void *pData, int
*pnWindowXSize, 
int *pnWindowYSize, int *pnPixelBuffer);
/**
* \brief Apply a raster processing function to each sub-window of a raster.
* 
* The input raster dataset is read in chunks of nWindowXSize * nWindowYSize
and each chunk is passed to the processing
* function. The output array from the function is written to the destination
dataset. 
* An optional 'pixel buffer' can be specified to allow overlaps between
successive windows. This is useful for 
* some algorithms, e.g. blob extraction, watershed/stream flow analysis,
convolution etc.
* Process specific data can be passed (e.g. configuration parameters). This
data is simply passed straight through to the processing
* function on each call. 
*
* @param processFn A GDALRasterProcessFn to apply to each sub window of the
raster. 
* 
* @param poSrcDataset The source raster dataset from which pixel values are
read
*
* @param poDstDataset The destination raster dataset to which pixel values
are written. Must support RasterIO in write mode.
*
* @param pData Process-specific data. This is passed straight through to the
GDALRasterProcessFn on each call. 
*
* @param pnWindowXSize The desired width of each read window. If NULL it
defaults to the 'natural' block size of the raster
*
* @param pnWindowYSize The desired height of each read window. If NULL it
defaults to the 'natural' block size.
*
* @param pnPixelBuffer A pixel buffer to apply to the read window. The read
window is expanded by pnPixelBuffer pixels in all directions such that 
*each window overlaps by pnPixelBuffer pixels. Ignored if NULL or 0
*
* @return a CPLErr enum indicating whether the function succesfully ran. 
*/


CPLErr GDALRunRasterProcesses(GDALRasterProcessFn *paProcessFn, int
nProcesses, GDALDataset *poSrcDataset, 
GDALDataset *poDstDataset, void **paData, int
*pnWindowXSize, 
int *pnWindowYSize, int *pnPixelBuffer);

/**
* \brief Apply multiple raster processing functions to each sub-window of a
raster
*
* For each window, the functions defined by the paProcessFn array are called
in turn, with the array output of the previous function forming the input
* to the next function. This allows processing 'toolchains' to be built
without having to create intermediate datasets, which can be less efficient
in time and space. 
*
*
* @param paProcessFn An array of GDALRasterProcessFn to apply to each sub
window of the raster
*
* @param nProcesses The size of paProcessFn
*
* @param poSrcDataset The source raster dataset from which pixel values are
read
*
* @param poDstDataset The destination raster dataset to which pixel values
are written. Must support RasterIO in write mode.
*
* @param paData an array of process-specific data objects of size
nProcesses. Each data object will be passed to the corresponding
GDALRasterProcessFn
*
* @param pnWindowXSize The desired width of each read window. If NULL it
defaults to the 'natural' block size of the raster
*
* @param pnWindowYSize The desired height of each read window. If NULL it

Re: [gdal-dev] GDAL raster processing library?

[jramm]

I can see where some similarities with other new and existing GDAL work could
be a blocker on this, but I also think this adds a a degree more flexibility
allowing potentially any kind of complex processing to be carried out
without worrying/bothering about boilerplate. 

It would be good to find some kind of middle ground where this flexibility
can be added in without duplicating other efforts?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/GDAL-raster-processing-library-tp5275948p5276078.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Error in GDALWarp to NWT_GRD

[jramm]

I have a strange error when reprojecting a dataset to NWT_GRD format. 

e.g.

gdalwarp -of NWT_GRD -ot Float32 -t_srs EPSG:4326 test.tif test.grd

very quickly returns


Stepping through, I've traced this to a call to `GetLockedBlockRef` in
IRasterIO (rasterio.cpp) which is called when GDALWarp attempts to write the
first warped 'block' to the GRD dataset. 

Somehow, a flag - `bJustInitialize` gets set to false, causing
GetLockedBlockRef to attemp a read operation before writing the block. As no
data has yet been written to the dataset, this read operation fails.
(Specifically, line 316 of grddataset.cpp in frmts/northwood)

Im struggling to find the cause of this. I am wondering whether I have
missed some implementation detail in the Create method of the NWT_GRD
driver?



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Error-in-GDALWarp-to-NWT-GRD-tp5276136.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Error in GDALWarp to NWT_GRD

[jramm]

jramm wrote
> gdalwarp -of NWT_GRD -ot Float32 -t_srs EPSG:4326 test.tif test.grd
> 
> very quickly returns

Sorry this meant to read:

very quickly returns:

0ERROR 1: /home/jamesramm/test.grd, band 1: IReadBlock failed at X offset 0,
Y offset 0
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Error-in-GDALWarp-to-NWT-GRD-tp5276136p5276137.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Error in GDALWarp to NWT_GRD

[jramm]

I added the following to the end of the Create method in
frmts/northwood/grddataset.cpp:


vsi_l_offset nFileSize = 1024 + nXSize * nYSize * 2;
if (VSIFTruncateL(poDS->fp, nFileSize) != 0) {
CPLError(CE_Failure, CPLE_FileIO,
"Failed to allocate space for GRD file");
delete poDS;
return NULL;
}
poDS->FlushCache(); // Write the header to disk.

Unfortunately still receiving the error. I wonder if it would be better if I
explicitly write the zeros with VSIFWriteL?




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Error-in-GDALWarp-to-NWT-GRD-tp5276136p5276347.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Memory increase in postgis raster driver

[jramm]

I am seeing a continuous, linear memory increase when reading sequential
windows of a large postgis dataset. This is the same regardless of my
GDAL_CACHEMAX settings
If I run he same code on an equivalent GTiff dataset, I do not see the same
increase in memory. 

I can see that the postgis dataset RasterIO implementation calls CacheTile
for each tile it reads, which I think uses the GDAL caching mechanism - is
it possible that, since each tile will be *new* (due to the sequential read
over the dataset), it is somehow cached outside of the gdal cache/ignoring
the CACHEMAX setting?

I'm just stabbing in the dark here as I can't really find a reason for it
yet



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Memory-increase-in-postgis-raster-driver-tp5277284.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Call for discussion on RFC 63 : Sparse datasets improvements

[jramm]

Forgive me a newbie question but I'm not so familiar with the RFC process:

- Is there any kind of timescale for this appearing in trunk or what needs
to happen for that to happen?




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-dev-Call-for-discussion-on-RFC-63-Sparse-datasets-improvements-tp5275443p5277676.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] GDAL raster processing library?

[jramm]

Hi Ari,
I began some work to clarify my ideas here: www.github.com/JamesRamm/GeoAlg
perhaps there is potential for merging with your project? For
neighbourhoods I provide 2 iterators - one simple block based and one
'buffered'. The block can be the natural block or user defined. I also
planned to provide a 'mosaic' iterator to support the planned 'reduce'
function. For callbacks I initially planned a simple function, but then
realised that many algorithms may wish to save state, so I changed to a
simple interface class where there is a single function to override. My
final plan was to expose opencv and other functionality (eg maybe taudem)
as a library of ready made callbacks.

If you think it worthwhile and possible, I would be open to porting
relevant parts to your code. Features I see as absolutely essential are:
- 'map' and 'reduce' functions for raster datasets, which handle large data
and user defined callbacks
- ability to save state in callbacks. This is essential for more complex
algorithms.
- A rich set of block iterators

A ready made library of a number of callbacks (preferably by exposing
functionality from existing well known libraries) would also be useful and
help build interest I think.

I will be in Bonn at foss4g from the Tuesday evening, it would be good to
meet!

On 8 Aug 2016 1:56 p.m., "Ari Jolma-2 [via OSGeo.org]" <
ml-node+s1560n5279906...@n6.nabble.com> wrote:

> Hi,
>
> Sorry for missing this last month due to holidays.
>
> In fact in my approach I do not want to limit to single pixel operations
> but support neighborhood operations too. That is necessary for watershed
> analysis for example. Making sure that the neighborhood is readily
> available for computations is probably what is causing the biggest problem
> for the code. I also have a callback support. See for example
> https://github.com/ajolma/gdal/blob/trunk/gdal/map_algebra/test.cpp where
> on the lines 5 to 12 is a callback function.
>
> When I left the code for my extended holidays I was doing profiling and
> intended to test another approach for caching. Now the caching is based on
> the native GDAL blocks, which is efficient for going through a single
> dataset. Since the block may be different for different datasets this
> approach may not be the best when there are more than one dataset. For
> example the raster algebra in QGIS uses an approach where the block is the
> same for all datasets. That simplifies the code and may be very efficient
> as a whole.
>
> I'll be at the code sprint in FOSS4G 2016 and that would be a good place
> to discuss the RFC too.
>
> Best,
>
> Ari
>
> 13.07.2016, 10:09, James Ramm kirjoitti:
>
> Peter,
> I think 'Grid Algebra' would be what Ari Jolma is proposing here:
> https://trac.osgeo.org/gdal/wiki/rfc62_raster_algebra
>
> As Even pointed out, there is some overlap, though my proposal is
> technically very different.
> The key differences I see are:
>
> - Users can submit functions which operate on each sub window of the
> raster, rather than an algebraic expression. This potentially allows for
> much more complicated algorithms to be used (e.g. I dont think it would be
> possible to run a watershed segmentation with a raster algebra
> implementation, or to have algorithms which behave differently depending on
> the 'location' within the raster or the values of surrounding pixels etc
> etc).
>
> - Functions can be chained together for a complex processing toolchain.
> Some overlap with VRT here, although again this introduces a little more
> flexibility.
>
> On 12 July 2016 at 15:47, Peter Halls <[hidden email]
> > wrote:
>
>> James,
>>
>>   in reality, are you not requesting an implementation of Tomlin's
>> 'Grid Algebra' in GDAL?  That defines the whole range of functions from
>> whole raster to pixel and has the distinct advantage of being both
>> published and extremely well known because of other implementations ...
>> which also provide ready-made reference bases for the GDAL implementors ...
>>
>> Best wishes,
>> Peter
>>
>> On 12 July 2016 at 15:39, James Ramm <[hidden email]
>> > wrote:
>>
>>> Hi Even,
>>>
>>> The difference I see with pixel functions is that, as far as I
>>> understand, the pixel function is applied per pixel, so there is no
>>> possibility of e.g. the pixel buffer when have the function apply to
>>> 'blocks'.
>>> I may be way off, but many of the algorithms we deal with require some
>>> kind of neighbourhood search - a polygonise algorithm or flow direction
>>> algorithm being good examples.
>>> I dont think VRT pixel functions allow this?
>>>
>>> So in that sense I'd see a VRT being 'just' another potential input data
>>> source.
>>>
>>> Perhaps VRT pixel functions could be expanded to also allow 'window'
>>> functions?
>>>
>>> A downside is it requires creating a VRT even when you only want to
>>> apply a such a function to a single dataset

Re: [gdal-dev] Memory increase in postgis raster driver

[jramm]

Hi 
We have been looking into this problem some more.

So it seems we can reproduce this by using postgis raster tiles of 100 x
100.
raster2pgsql cmd was:
raster2pgsql bigtiff.tif -C -r -s 27700 -t 100x100 -P -I -M -Y
 | psql 

The effect seems more pronounced using GDAL 2.2 (trunk) vs GDAL 2.1

With a tilesize of 256*256 it is fine. 

This is a very large tif btw; ~15 x 20 pixels





--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Memory-increase-in-postgis-raster-driver-tp5277284p5280563.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] VRT derived band pixel functions written in Python

[jramm]

Trying to run this using a function relying on scipy.ndimage...

When running gdal_translate on the VRT, I get ImportError: No module named
scipy.ndimage
This comes after successfully import numpy. scipy.ndimage will happily
import within the python interpreter.
Any tips on how to track this down/debug? 

The entire VRT file is as follows:


  PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB
1936",DATUM["OSGB_1936",SPHEROID["Airy
1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",40],PARAMETER["false_northing",-10],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","27700"]]
  100180.0,5.0,0.0,1215730.0,0.0,-5.0
  

  F:\tif_data\large_sparse.tif
  
  
4
TRUE
  

extract_blobs
Python


5

  





--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-dev-VRT-derived-band-pixel-functions-written-in-Python-tp5285323p5285882.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] VRT derived band pixel functions written in Python

[jramm]

Trying to run this using a function relying on scipy.ndimage...

When running gdal_translate on the VRT, I get ImportError: No module named
scipy.ndimage
This comes after successfully import numpy. scipy.ndimage will happily
import within the python interpreter.
Any tips on how to track this down/debug? 

The entire VRT file is as follows:


  PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB
1936",DATUM["OSGB_1936",SPHEROID["Airy
1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],AUTHORITY["EPSG","6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4277"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",40],PARAMETER["false_northing",-10],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","27700"]]
  100180.0,5.0,0.0,1215730.0,0.0,-5.0
  

  F:\tif_data\large_sparse.tif
  
  
4
TRUE
  

extract_blobs
Python


5

  





--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-dev-VRT-derived-band-pixel-functions-written-in-Python-tp5285323p5285883.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev