RE: [gdal-dev] netCDF to Shapefile
Ok, here are more details: I am using Python (2.6) together with netCDF4, numpy and gdal/ogr. So far, I am able to read the dimensions, attributes and variables. I have to choose records depending on certain time slices and export these records together with their dimensions into a shapefile. So, if four time slices are selected four shp files will be created. Currently, I am designing a GUI with Tkinter for data I/O and time slice selection. For the conversion task i would prefer a direct method without caching into a file. I may use normal arrays or pytables Just wanted to know what is the more efficient method for caching the data. Thanks for your support Mike! Date: Wed, 16 Mar 2011 16:21:23 +1100 Subject: Re: [gdal-dev] netCDF to Shapefile From: mdsum...@gmail.com To: manuelrai...@hotmail.com CC: gdal-dev@lists.osgeo.org Hi Manuel, I think you will still need to be more specific. You've indicated that the data is points, but not what the variables in the file are, which variables need to be interrogated for coordinates vs. values, what their dimensions are, whether you want several shapefiles, or just one with (say) multiple attributes, etc. There are a number of ways to convert such data to shapefile, but it might be best to choose a reasonably easy scripting language to do it, such as Python or R - and so there are probably better forums than this for your problem. Either way you'll need to provide a lot more detail - it's likely to be a two step process, with converting the NetCDF grids to tables (in memory or file), and then those to point shapefiles. Cheers, Mike. On Wed, Mar 16, 2011 at 3:17 PM, Manuel Rainer manuelrai...@hotmail.com wrote: Thanks for your help, it has to be a Shapefile because we do some further processing (symbology, map creation...). The netCDF file contains hydrodynamics point data. Best regards Manuel Date: Wed, 16 Mar 2011 09:32:16 +0530 Subject: Re: [gdal-dev] netCDF to Shapefile From: chaitanya...@gmail.com To: manuelrai...@hotmail.com CC: gdal-dev@lists.osgeo.org Manuel, Shapefile is for vector formats alone. Try a raster format to convert the netCDF data. As for the performance, you don't need to cache anything. But it may help to tweak the GDAL_CACHEMAX setting. http://trac.osgeo.org/gdal/wiki/FAQRaster#Howtoimprovegdalwarpperfomance On Wed, Mar 16, 2011 at 8:21 AM, Manuel Rainer manuelrai...@hotmail.com wrote: Hi, I would like to convert certain time slices (records) out of a netCDF file to a Shapefile. The netCDF file has 7 GB, so what is the best proceed to get an appropriate performance? I mean, should the data be cached (e.g. in an ASCII file, ...) before writing a Shapefile or can the shp directly be generated? Perhaps you have some experiences... Thanks for your help! Best regards Manuel ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Best regards, Chaitanya kumar CH. /tʃaɪθənjə/ /kʊmɑr/ +91-9494447584 17.2416N 80.1426E ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Michael Sumner Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia e-mail: mdsum...@gmail.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] netCDF to Shapefile
On Wed, Mar 16, 2011 at 3:51 AM, Manuel Rainer manuelrai...@hotmail.com wrote: Hi, I would like to convert certain time slices (records) out of a netCDF file to a Shapefile. The netCDF file has 7 GB, so what is the best proceed to get an appropriate performance? I mean, should the data be cached (e.g. in an ASCII file, ...) before writing a Shapefile or can the shp directly be generated? Perhaps you have some experiences... Hi I don't know what you mean with performance (maybe the time needed to make the export or the subsequent access to them?), but as far as I remember shapefile has a file size limit of 2 Gb, I personally would avoid that approach. Here PostGis seems to well suit in the scenario, but if you don't have the possibility to access/create a PostGis infrastructure, I personally would give a try to Spatialite. If performance are really a concern, you could instead opt for a NoSQL database: CouchDb has a nice spatial module [0] as far as I have heard. You could use the shp2geocouch approach utility [1] (first export to json with ogr2ogr then json to couchdb) for loading your data. Please note that if you opt for the latest, I am not sure if you will have tools for symbology and map creation at this time, neither you can use the GDAL API. best regards Paolo [0] http://vmx.cx/blog/2009-11-17/geodata-and-couchdb.htm [1] https://github.com/maxogden/shp2geocouch/blob/master/bin/shp2geocouch -- Paolo Corti Geospatial software developer web: http://www.paolocorti.net twitter: @paolo_corti ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Calculate footprints of shapefiles
Chaitanya thanks for the hints. So a sort of Union operation seems to be necessary. I need to check how fast that is with OGR on larger shapefiles (I want to store in PostGIS just the footprint, not the full geometry collection). The ConvexHull could be a good alternative to a simplified boundary. I understand how the normal Union() works, it merges 2 geometries into a new one. But I don't get it what UnionCascade() is doing. It applies a Union on a single geometry and returns the new geometry, what is it actually unioning/merging then? Regards, Armin Original-Nachricht Datum: Wed, 16 Mar 2011 10:06:57 +0530 Von: Chaitanya kumar CH chaitanya...@gmail.com An: armin.bur...@gmx.net CC: Marius Jigmond mariusjigm...@hotmail.com, gdal-dev@lists.osgeo.org Betreff: Re: [gdal-dev] Calculate footprints of shapefiles Armin, This can be done in either OGR or PostGIS. For OGR, refer to the OGRGeometry class reference, especially the ConvexHull(), Union() and UnionCascaded() functions. For PostGIS, ST_Boundary() and ST_Union(). -- GMX DSL Doppel-Flat ab 19,99 Euro/mtl.! Jetzt mit gratis Handy-Flat! http://portal.gmx.net/de/go/dsl ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] netCDF to Shapefile
I have done this but not to a shapefile direct, but into a postgis database. Use for loops to loop through the temporal dimension, then the lon and lat as sub loops. Then you create an insert string to input each line of data into the sql database. In our NetCDF we have 720x360 observations spatially and 1308 temporal dimension observations. Thats a lot of data! Andreas 2011/3/16 Paolo Corti pco...@gmail.com On Wed, Mar 16, 2011 at 3:51 AM, Manuel Rainer manuelrai...@hotmail.com wrote: Hi, I would like to convert certain time slices (records) out of a netCDF file to a Shapefile. The netCDF file has 7 GB, so what is the best proceed to get an appropriate performance? I mean, should the data be cached (e.g. in an ASCII file, ...) before writing a Shapefile or can the shp directly be generated? Perhaps you have some experiences... Hi I don't know what you mean with performance (maybe the time needed to make the export or the subsequent access to them?), but as far as I remember shapefile has a file size limit of 2 Gb, I personally would avoid that approach. Here PostGis seems to well suit in the scenario, but if you don't have the possibility to access/create a PostGis infrastructure, I personally would give a try to Spatialite. If performance are really a concern, you could instead opt for a NoSQL database: CouchDb has a nice spatial module [0] as far as I have heard. You could use the shp2geocouch approach utility [1] (first export to json with ogr2ogr then json to couchdb) for loading your data. Please note that if you opt for the latest, I am not sure if you will have tools for symbology and map creation at this time, neither you can use the GDAL API. best regards Paolo [0] http://vmx.cx/blog/2009-11-17/geodata-and-couchdb.htm [1] https://github.com/maxogden/shp2geocouch/blob/master/bin/shp2geocouch -- Paolo Corti Geospatial software developer web: http://www.paolocorti.net twitter: @paolo_corti ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Pickling OGR Geometry
Hi all, I'm looking at pickling ogr geometry. I took the following code from the gdal autotests to try it out: geom_wkt = 'MULTIPOLYGON( ((0 0,1 1,1 0,0 0)),((0 0,10 0, 10 10, 0 10),(1 1,1 2,2 2,2 1)) )' geom = ogr.CreateGeometryFromWkt(geom_wkt) p = pickle.dumps(geom) g = pickle.loads(p) if not geom.Equal(g): print pickled geometries were not equal else: print pickled geometries equal The geometries are equal and appears to work fine but I get an error printed out right after the loads call: ERROR 1: Empty geometries cannot be constructed If I do ogr.UseExceptions() I get the stack trace: Traceback (most recent call last): File test.py, line 32, in module g = pickle.loads(p) File /usr/lib/python2.6/pickle.py, line 1374, in loads return Unpickler(file).load() File /usr/lib/python2.6/pickle.py, line 858, in load dispatch[key](self) File /usr/lib/python2.6/pickle.py, line 1133, in load_reduce value = func(*args) File /usr/local/lib/python2.6/dist-packages/GDAL-1.7.2-py2.6-linux-x86_64.egg/osgeo/ogr.py, line 2935, in __init__ this = _ogr.new_Geometry(*args, **kwargs) RuntimeError: Empty geometries cannot be constructed Is this a problem or is it just noise? I'm using GDAL 1.7.2. Thanks, Jason ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Need help generating output with alpha channels using Python API
On 11-03-14 12:50 PM, Michal Migurski wrote: Hello, I'm having some difficulties understanding how to get GDAL to generate images with usable alpha channels. I have 3-channel RGB input JPEG image, delivered to GDAL as a VRT with a projection and ground control points, which I'm warping to an output that's no longer rectangular, leaving background areas exposed. Here is an example output: http://things.teczno.com/gnomotile.png The black parts are intended to be transparent, but I haven't been able to understand how to make that work. Here's the relevant part of my Python code: http://dpaste.com/hold/500308/ I've tried to switch to a four channel output which gets me what I think are CMYK channels. I've tried to use SetNoDataValue on the destination bands to make the background purple so it can be easily knocked out. I've tried to create the GTiff output dataset using the ALPHA=YES creation option, but it seemingly doesn't make a difference. None of these ideas has worked - does anyone have any ideas on how the Python API can be used to create transparent output? Mike, I would have thought preinitializing the destination dataset to a marker value should have done the trick. I'm not sure what you mean by using SetNoDatavalue(). Are you aware that this just records metadata with the band indicating what the nodata value should be? Normally I would suggest adding an alpha band on the output and ensuring it is used. But I see GDALReprojectImage() does not make it easy to do this particularly the way it is bound in Python which does not allow one to pass in the warp options structure. So one approach is for you to precreate the output image, initializing with some particular nodata pixel value and then using that in post processing. Another approach might be for the me to modify GDALReprojectImage to utilize an alpha band properly if it exists on the destination image. I can do this in trunk if that would be helpful to you. Best regards, -- ---+-- I set the clouds in motion - turn up | Frank Warmerdam, warmer...@pobox.com light and sound - activate the windows | http://pobox.com/~warmerdam and watch the world go round - Rush| Geospatial Programmer for Rent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_translate format listing is too long
Hello all, Question - is there some reason why gdal_translate needs to show a format listing when invoked on the command line with no arguments? On my system this outputs 61 lines which is nearly a maximized shell on my 24 monitor. Wouldn't requiring the --formats flag be more consistent with other tools? I can file a feature request for this - just though I'd ask first. Cheers, Jamie ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Get pixel values from all bands
Hi all
I work on raster processing application (developed with C++ and Qt).
For reading rasters we use GDAL and it's C++ API.
For one task we need to get pixel values from all raster bands. For
this I get all GDALRasterBands and store them in list. Then in loop
read image row by row and store row from each band in another list.
Here is simplified code
// iterate over image row by row
for ( int row = 0; row ySize; ++row )
{
// read one row from each band and store it in list
for ( int b = 0; b bandCount; ++b )
{
float *scanline;
scanline = (float *) CPLMalloc( sizeof( float ) * xSize );
bands[ i ]-RasterIO( GF_Read, 0, row, xSize, 1, scanline, xSize,
1, GDT_Float32, 0, 0 );
rows[ i ] = scanline;
}
// iterate on pixels in row
for (int col = 0; col ySize; ++col)
{
// get pixel value from bands
float pixel[bandCount];
for ( int r = 0; r bandCount; ++r )
{
// get value in (row, col) from band r
double pixel[r] = rows[ r ][ col ];
}
// do something with this pixel
}
}
Is this correct approach? Or maybe there is another more correct way
to do this?
Also I'm not sure where I should free memory allocated with CPLMalloc().
After reading pixels and before reading next rows or at end when
all rows processed?
I'm not a big C++ programmer, so any help is very appreciated
Thanks
--
Alexander Bruy
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Get pixel values from all bands
Alexander Bruy wrote:
Here is simplified code
// iterate over image row by row
for ( int row = 0; row ySize; ++row )
{
// read one row from each band and store it in list
for ( int b = 0; b bandCount; ++b )
{
float *scanline;
scanline = (float *) CPLMalloc( sizeof( float ) * xSize );
bands[ i ]-RasterIO( GF_Read, 0, row, xSize, 1, scanline, xSize,
1, GDT_Float32, 0, 0 );
rows[ i ] = scanline;
}
// iterate on pixels in row
for (int col = 0; col ySize; ++col)
{
// get pixel value from bands
float pixel[bandCount];
for ( int r = 0; r bandCount; ++r )
{
// get value in (row, col) from band r
double pixel[r] = rows[ r ][ col ];
}
// do something with this pixel
}
}
Is this correct approach? Or maybe there is another more correct way
to do this?
Seems OK to me, except that I don't understand what is the cache of row
array for. Why not to process pixels straight from scanline.
Also I'm not sure where I should free memory allocated with CPLMalloc().
After reading pixels and before reading next rows or at end when
all rows processed?
Don't bother with CPLMalloc, let the compiler to clean after you.
typedef std::vectorfloat raster_data_t;
raster_data_t scanline(xSize);
bands[ i ]-RasterIO( GF_Read, 0, row, xSize, 1, scanline[0], xSize,
1, GDT_Float32, 0, 0 );
Best regards,
--
Mateusz Loskot
http://mateusz.loskot.net
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Get pixel values from all bands
Alexander,
Not sure if it meets you situation, but there is a command line utility
with similar functionality, gdallocationinfo,
http://gdal.org/gdallocationinfo.html
HTH, Eli
On 3/16/2011 at 1:21 PM, in message
20110316222128.b943f4e3.alexander.b...@gmail.com, Alexander Bruy
alexander.b...@gmail.com wrote:
Hi all
I work on raster processing application (developed with C++ and Qt).
For reading rasters we use GDAL and it's C++ API.
For one task we need to get pixel values from all raster bands. For
this I get all GDALRasterBands and store them in list. Then in loop
read image row by row and store row from each band in another list.
Here is simplified code
// iterate over image row by row
for ( int row = 0; row ySize; ++row )
{
// read one row from each band and store it in list
for ( int b = 0; b bandCount; ++b )
{
float *scanline;
scanline = (float *) CPLMalloc( sizeof( float ) * xSize );
bands[ i ]-RasterIO( GF_Read, 0, row, xSize, 1, scanline, xSize,
1, GDT_Float32, 0, 0 );
rows[ i ] = scanline;
}
// iterate on pixels in row
for (int col = 0; col ySize; ++col)
{
// get pixel value from bands
float pixel[bandCount];
for ( int r = 0; r bandCount; ++r )
{
// get value in (row, col) from band r
double pixel[r] = rows[ r ][ col ];
}
// do something with this pixel
}
}
Is this correct approach? Or maybe there is another more correct way
to do this?
Also I'm not sure where I should free memory allocated with CPLMalloc().
After reading pixels and before reading next rows or at end when
all rows processed?
I'm not a big C++ programmer, so any help is very appreciated
Thanks
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Reproject Globcover
Greetings I pretend to reproject a known JRC product called GLOBCOVER to UTM WGS84 (zone for France) using gdalwarp keeping spatial resolution (300m) and null_values. Can anyone give me a suggestion on how to reproject it using gdal? Bellow is the gdalinfo information Driver: GTiff/GeoTIFF Files: d:\data\Local\GLobcover\1.tif Size is 18720, 17640 Coordinate System is: GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.2572235630016, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4326]] Origin = (-11.00137062376,83.00138883284) Pixel Size = (0.0027808,-0.0027780) Metadata: AREA_OR_POINT=Area Image Structure Metadata: COMPRESSION=LZW INTERLEAVE=BAND Corner Coordinates: Upper Left ( -11.0013889, 83.0013889) ( 11d 0'5.00W, 83d 0'5.00N) Lower Left ( -11.0013889, 34.0013889) ( 11d 0'5.00W, 34d 0'5.00N) Upper Right ( 40.9986111, 83.0013889) ( 40d59'55.00E, 83d 0'5.00N) Lower Right ( 40.9986111, 34.0013889) ( 40d59'55.00E, 34d 0'5.00N) Center ( 14.9986111, 58.5013889) ( 14d59'55.00E, 58d30'5.00N) Band 1 Block=18720x1 Type=Byte, ColorInterp=Gray NoData Value=0 Metadata: BAND=CL6 Thanks ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Reproject Globcover
On Wed, Mar 16, 2011 at 5:12 PM, Monica Buescu monicabuescu1...@gmail.com wrote: I pretend to reproject a known JRC product called GLOBCOVER to UTM WGS84 (zone for France) using gdalwarp keeping spatial resolution (300m) and null_values. Can anyone give me a suggestion on how to reproject it using gdal? Try: gdalwarp -t_srs EPSG:32631 -tr 300 300 -dstnodata 0 input file output file EPSG 32631 is WGS84 UTM zone 31N. See http://spatialreference.org/ref/epsg/32631/ -b -- Brian Claywell bcclayw...@gmail.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Calculate footprints of shapefiles
thanks for the hints. So a sort of Union operation seems to be necessary. I need to check how fast that is with OGR on larger shapefiles (I want to store in PostGIS just the footprint, not the full geometry collection). The Convex Hull could be a good alternative to a simplified boundary. it depends on your shapes -- if you have a very concave shape, a convex hull can be a poor representation. You might try simplifying before the Union - that would make the Union faster. -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/ORR(206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception chris.bar...@noaa.gov ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: Re: [gdal-dev] Where is this nmake.opt file located?
Natasha, After you install the GDAL and GDAL-Oracle packages, you can access the command line applications in the OSGeo4W Shell command window. http://trac.osgeo.org/osgeo4w/#QuickStartforOSGeo4WUsers You can find the info on connecting to Oracle Spatial in the GDAL website. For vector data access: http://www.gdal.org/ogr/drv_oci.html For raster data access: http://www.gdal.org/frmt_georaster.html On Thu, Mar 17, 2011 at 12:51 AM, natasha chatterjee natashachatter...@rediffmail.com wrote: Thanks I installed Osgeo4w but now I am not getting how to connect GDAL with Oracle spatial.Please tell me the process and what are the prerequisites to bechecked.(How to proceed?) Kindly help me. Regards Natasha Tynaks On Tue, 15 Mar 2011 23:57:06 +0530 wrote Natasha, nmake.opt file is for when you compile the source. It is in the root directory of the GDAL source tree. If you intend to use Oracle Spatial with GDAL on Windows, you are better off with OSGeo4W than FWTools. It provides oracle support as a separate plugin. http://trac.osgeo.org/osgeo4w/ On Tue, Mar 15, 2011 at 10:52 PM, natasha chatterjee wrote: Hi All I am using FWTools.Can anyone tell me where is this nmake.opt file located on our system?if its not there then how to get it?(As i want to configure oracle spatial with GDAL) please reply soon. Thanks and regards Natasha ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Best regards, Chaitanya kumar CH. /tʃaɪθənjə/ /kʊmɑr/ +91-9494447584 17.2416N 80.1426E http://sigads.rediff.com/RealMedia/ads/click_nx.ads/www.rediffmail.com/signatureline.htm@Middle? -- Best regards, Chaitanya kumar CH. /tʃaɪθənjə/ /kʊmɑr/ +91-9494447584 17.2416N 80.1426E ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Calculate footprints of shapefiles
Armin, Cascaded Union works pretty much like Union except that it is optimized to work on more than two geometries. Both OGR and PostGIS uses the GEOS library to perform this. Read this blog entry by Paul Ramsey: http://blog.cleverelephant.ca/2009/01/must-faster-unions-in-postgis-14.html On Wed, Mar 16, 2011 at 3:41 PM, Armin Burger armin.bur...@gmx.net wrote: Chaitanya thanks for the hints. So a sort of Union operation seems to be necessary. I need to check how fast that is with OGR on larger shapefiles (I want to store in PostGIS just the footprint, not the full geometry collection). The ConvexHull could be a good alternative to a simplified boundary. I understand how the normal Union() works, it merges 2 geometries into a new one. But I don't get it what UnionCascade() is doing. It applies a Union on a single geometry and returns the new geometry, what is it actually unioning/merging then? Regards, Armin Original-Nachricht Datum: Wed, 16 Mar 2011 10:06:57 +0530 Von: Chaitanya kumar CH chaitanya...@gmail.com An: armin.bur...@gmx.net CC: Marius Jigmond mariusjigm...@hotmail.com, gdal-dev@lists.osgeo.org Betreff: Re: [gdal-dev] Calculate footprints of shapefiles Armin, This can be done in either OGR or PostGIS. For OGR, refer to the OGRGeometry class reference, especially the ConvexHull(), Union() and UnionCascaded() functions. For PostGIS, ST_Boundary() and ST_Union(). -- GMX DSL Doppel-Flat ab 19,99 Euro/mtl.! Jetzt mit gratis Handy-Flat! http://portal.gmx.net/de/go/dsl -- Best regards, Chaitanya kumar CH. /tʃaɪθənjə/ /kʊmɑr/ +91-9494447584 17.2416N 80.1426E ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
