Re: [gdal-dev] Processing builder
Hi Frank, your question is not about GDAL but about QGIS and should be asked in the QGIS-users mailing list. Your issue is that you are using online DEM source. You should first download it to the local machine and only then call Processing algorithm against saved file. Hope this helps. чт, 5 бер. 2020 о 10:26 Frank Miller пише: > > Hi, > I am very new in this group and have a question you may be very familiar > with. I am trying to add 'Clip raster by mask layer' tool to processing > modeler and download DEM data from an online resource but so far there is no > luck. > Does GDAL works in Processing Modeler? If yes what is my problem? > > Regards, > Frank > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/gdal-dev -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] About CMake build again
Hi Dmitry, thanks for explanation. I see that there is an option to use system libraries. Maybe I'm wrong, but this is also possible with current build system and CMake scripts from #7080. There is an option to build deps by myself and use them to build GDAL, but again this is also possible with current build system and CMake scripts from #7080. I agree that automatic fetching of missed dependencies is nice feature. And can simplify live in some cases. But there are questions: 1. what if dependency hosted on BitBucket, SourceForge or any self-hosted VCS and not in your repository? 2. what if dependency use build system different from CMake? 3. what if upstream does not want/does not see abny benefits in moving to GitHub and/or switching to CMake? Correct me if I'm wrong, but most of the borsh's benefits require changes not only in GDAL, but also in all upstream projects, as well as infrastructure changes for them. Or why you maintain forks for every lib needed by GDAL in your repository? If everything is fine and flexible why not use upstream projects directly? 2017-10-29 15:39 GMT+02:00 Dmitry Baryshnikov <bishop@gmail.com>: > 1. Build gdal with all dependencies getting them from github > (WITH_EXTERNAL). Preferable for Windows > > 2. Build gdal using the system libraries. Preferable for Linux > > 3. Build gdal using the dependency libraries build by user (out of source) > and registered in CMake package registry. This is I use now. This script > help me to clone all libraries, build them and register them in CMake > package registry > (https://github.com/nextgis-borsch/borsch/blob/master/opt/tools.py#L134). > > My opinion is that different options is much flexible and suits for many > scenarios as very limited solution from > https://trac.osgeo.org/gdal/ticket/7080. -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] About CMake build again
host system. > > One can see this is all very flexible. > > > What about GDAL. > > 1. After unification GDALDataset and OGRDatasource current sources tree is > not fit for this new logic of GDAL classes. I rearranged sources more closer > to GDAL classes and CMake needs. Main changes are moving raster and vector > drivers inside drivers folder > (https://github.com/nextgis-borsch/lib_gdal/tree/master/drivers).This > simplify situation where different drivers need the same dependency library > (libpg, libsqlite, etc.). Also there are several raster/vector drivers which > need a separate directory but now presented in ogr or frmts directories. > There are some bad decisions I made - for example I moved unit tests into > separate repository - this was a mistake. We will return unit tests back to > GDAL repository. > > An example of cmake friendly way see > https://github.com/nextgis-borsch/lib_gdal/blob/master/drivers/vector/CMakeLists.txt. > The driver developer must only create new folder and put CMakeLists.txt file > into it. The upper CMake script will find new driver and add it to GDAL > build. In common cases no need to modify upper CMake scripts. > > 2. I remove third-party code from drivers folders (TIFF, GeoTIF, PNG, JPEG > etc). All this code are in separate repositories. I don't see the difference > to get this code from git pull from main GDAL repository or from the > separate repository via find_anyproject process. Current GDAL repository > looks like the https://github.com/nextgis-borsch packed in one repository. > > > In conclusion: > > 1. Borsch added flexible and useful features and not remove existing > approach > > 2. The current cmaked GDAL are in production in my company more than a year > on Windows, Linix, Mac OS, iOS, Android. > > 3. I'm ready to discuss and improve current solution. Any help are welcome > > 4. I also will be happy to contribute directly or via PR into GDAL trunk > instead of backporting in both directions improvements that we do in GDAL . > > > Finally: > > Find the link to page with the CMake in GDAL discussion - > https://trac.osgeo.org/gdal/wiki/CMake > > -- > Best regards, > Dmitry > > > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/gdal-dev -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Adding more stuff to Python bindings
Hi devs, I want to add to Python bindings some functions which currently are available only in C API. Are there any instructions about how to proceed? As I understand, I need to add function declaration to the swing/include/Operations.i and also add missed type definitions somewhere. Is this correct approach? Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Crop raster without changing output pixel size
Hi all, I have a georeferenced raster and want to crop it using shapefile as mask layer. I use next command: gdalwarp -q -cutline myshape.shp -crop_to_cutline -of GTiff input.tif output.tif All fine except one thing: the output raster has another pixel size than input. For example, input raster has pixel size 0.0083,-0.0083, and output has 0.00833129,-0.00833387. I suspect this is related to rounding and 0.5 pixel offset (top left corner vs. pixel center). When mask is rectangle, I can use gdal_traslate to clip, but often I need to crop raster by non-rectangular mask. Is it possible to do this without raster grid resampling? Keeping pixel size is very important for further analysis. Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Crop raster without changing output pixel size
2013/9/25 Jukka Rahkonen jukka.rahko...@mmmtike.fi: Have a try by -tr together with -tap http://www.gdal.org/gdalwarp.html I guess it will give about what you want. Thanks for hints, I'll try this. But seems to use -tr option I should first calculate resolution of the output raster. This can work for time from time usage, but not for automated processing. Anyway, thanks for your help. -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] zonal statistics with gdal
Hi, you can use GDAL Python bindigs for this [0, 1]. The process may look as open shape and raster, then convert pixel coordinates into map (shape) coordinates and process only pixels inside your polygon. [0] http://gdal.org/gdal_tutorial.html [1] http://gdal.org/ogr/ogr_apitut.html 2012/6/19 jdmorgan jdmor...@unca.edu: Hello, Is there a way to do something along the lines of zonal statistics with gdal. Basically, I have a county shape file and a classified raster dataset. I would like to count the number of raster pixels/county (FIPS). Is this possible with GDAL/OGR? If so can someone point me in the right direction? If so, awesome as I plan to do some automation with python. Thanks, Derek Hope this helps -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Request to Visual Studio 2008+ users
Hi, just finished compiling using your instructions and it compiles fine 2012/4/22 Mateusz Loskot mate...@loskot.net: Folks, I'm looking for a but of assistance from GDAL hackers who use Visual Studio 2008 or later. I'd be thankful if someone could conduct the following test for me against the current GDAL trunk: 0) Grab sources from GDAL from SVN trunk or *clean* the sources tree thoroughly, so cpl_config.h as well as all intermediate, .obj and output binaries are wiped out 1) Create nmake.local 2) Put the following 6 lines in nmake.local !IFNDEF DEBUG OPTFLAGS= $(CXX_ANALYZE_FLAGS) /nologo /MD /EHsc /Ox /D_CRT_SECURE_NO_DEPRECATE /D_CRT_NONSTDC_NO_DEPRECATE /DNDEBUG !ELSE OPTFLAGS= $(CXX_ANALYZE_FLAGS) /nologo /MDd /EHsc /Zi /D_CRT_SECURE_NO_DEPRECATE /D_CRT_NONSTDC_NO_DEPRECATE /Fd$(GDAL_ROOT)\gdal$(VERSION).pdb /DDEBUG !ENDIF 3) nmake /f makefile.vc MSVC_VER= DEBUG-1 where is 1500 (VS2008) or 1600 (VS2010). 4) Does GDAL build successfully for you? If not, what errors are reported. Thanks in advance! Best regards, -- Mateusz Loskot, http://mateusz.loskot.net -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Request to Visual Studio 2008+ users
2012/4/22 Mateusz Loskot mate...@loskot.net: On 22 April 2012 14:24, Alexander Bruy alexander.b...@gmail.com wrote: Hi, just finished compiling using your instructions and it compiles fine Would you mine sharing details of Visual Studio version you use? I use Microsoft Visual Studio 2008 Version 9.0.30729.1 SP (Express Edition). Also I have Windows Server 2003 R2 Platform SDK installed. -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal2tiles crashes on some rasters
Hi all, I have troubles using gdal2tiles. When I try to generate tiles from raster it fails with error Generating Base Tiles: ERROR 5: Illegal values for buffer size ERROR 5: Illegal values for buffer size Traceback (most recent call last): File /usr/bin/gdal2tiles.py, line 2241, in module gdal2tiles.process() File /usr/bin/gdal2tiles.py, line 478, in process self.generate_base_tiles() File /usr/bin/gdal2tiles.py, line 1276, in generate_base_tiles dsquery.WriteRaster(wx, wy, wxsize, wysize, data, band_list=list(range(1,self.dataBandsCount+1))) File /usr/lib/python2.6/site-packages/GDAL-1.8.1-py2.6-linux-i686.egg/osgeo/gdal.py, line 746, in WriteRaster buf_pixel_space, buf_line_space, buf_band_space ) TypeError: not a string Here is my command gdal2tiles.py -z '8-14' --srcnodata=0 merc_SHL_tmv.tif ./output gdalinfo reports that raster georeferenced Coordinate System is: PROJCS[WGS 84 / UTM zone 54N, GEOGCS[WGS 84, DATUM[WGS_1984, SPHEROID[WGS 84,6378137,298.257223563, AUTHORITY[EPSG,7030]], AUTHORITY[EPSG,6326]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433], AUTHORITY[EPSG,4326]], PROJECTION[Transverse_Mercator], PARAMETER[latitude_of_origin,0], PARAMETER[central_meridian,141], PARAMETER[scale_factor,0.9996], PARAMETER[false_easting,50], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AUTHORITY[EPSG,32654]] When I set input CRS in gdal2tiles command and try again same error appears. Also I try to warp input raster to Google Mercator using this command gdalwarp -t_srs '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs' SHL_tmv.tif merc_SHL_tmv.tif and generate tiles from reprojected image. Same result. But for some other raster after reprojecting with same command gdal2tiles works. Any ideas what is wrong? Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Strange behavior of gdalwarp with -overwrite option
Hi all, in gdalwarp documentation [0] listed option -overwrite with next description: -overwrite: (GDAL = 1.8.0) Overwrite the target dataset if it already exists. But when I try to reproject image using next command gdalwarp -t_srs EPSG:4326 -overwrite mosaic.tif mosaic.tif i get next error Creating output file that is 9122P x 2341L. ERROR 4: `mosaic.tif' not recognised as a supported file format. and mosaic.tif becomes 8 bytes size file (seems only tif header). Is this bug or I misunderstand -overwrite option usage? Tested on Slackware 13.37, GDAL 1.8.1, released 2011/07/09. [0] http://www.gdal.org/gdalwarp.html Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Strange behavior of gdalwarp with -overwrite option
Hi, 2011/11/8 Even Rouault even.roua...@mines-paris.org: And you spend half and hour understanding why it has not been reprojected in EPSG:. The default behaviour being that once the target dataset exists, it is only updated, but not modified in depth. gdalwarp something.tif something.tif -t_srs EPSG: makes no sense, with or without -overwrite (and -overwrite will actually kill your target, ahem, source - since source=target - as you have seen !). Reprojection in place cannot work. You must specify a different target dataset. I can imagine that a safety check would make sense to detect such obvious error... Thanks for clarification. So answer on my question is: in place reprojection not working and -overwrite option only necessary if I want to recreate output image (not usefull IMHO, because I can specify another one output). -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] GDAL Pascal bindings
Hi all, some time ago I start to write Pascal bindings for GDAL library. Here is current results, maybe anyone find them useful http://gis-lab.info/share/alexbruy/gdal_pascal.tar.bz2 Currently only OGR part tested (see examples in demo directory). I continue to work on it, but there are also several projects I'm involved, so any help is welcome. Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Raster size and ReadAsArray()
Hi, There is a well-know problem: reading really large rasters or bands into memory with DataSource.ReadAsArray() method impossible due memory limitations. For example, when I try to read one band with size 53109x29049 I get error: File C:\OSGeo4W\apps\Python25\lib\site-packages\osgeo\gdal.py, line 727, in ReadAsArray return gdalnumeric.DatasetReadAsArray( self, xoff, yoff, xsize, ysize, buf_obj ) File C:\OSGeo4W\apps\Python25\lib\site-packages\osgeo\gdal_array.py, line 162, in DatasetReadAsArray return BandReadAsArray( ds.GetRasterBand(1), xoff, yoff, xsize, ysize, buf_obj = buf_obj) File C:\OSGeo4W\apps\Python25\lib\site-packages\osgeo\gdal_array.py, line 228, in BandReadAsArray ar = numpy.empty([buf_ysize,buf_xsize], dtype = typecode) ValueError: dimensions too large. I want to know is it possible to get maximum raster size that can be handled using ReadAsArray() without errors because I want to implement a fallback algorithm for large rasters in my tool. Currently I try to implement it with try-except statement but maybe there is more elegant solution? Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Raster size and ReadAsArray()
Hi 2011/8/3 Antonio Valentino antonio.valent...@tiscali.it: In my experience using too large chunks of memory can cause, paradoxically, slowdowns. My suggestion is to define a reasonable maximum size for arrays in your application/library and switch to the fallback algorithm for large rasters every time that MAX_SIZE is exceeded, even if using ReadAsArray still works. Thanks for suggestion, I already think about this. 2011/8/3 Even Rouault even.roua...@mines-paris.org: It is difficult to predict in all cases. But there are a few things to keep in mind. Do you use a 32 bit build of GDAL or a 64 bit one ? If it is a 32 bit, then your memory allocations are generally limited to 2 GB by process, and much of the time even less because it is difficult to get reliably a continuous area of virtual memory of such a size. 53109 x 29049 = 1.5 GB if your data is of type Byte and if your dataset has only one band (multiply by the number of bands and the size in byte of the band data type). With a 64 bit build of GDAL and enough RAM, that should work more reliably. Personally I use 32-bit system and 32-bit GDAL, but I want to make my tool flexible as much as possible. This is why I look for method to get maximum raster size that fits into memory limits. Currently I try to implement it with try-except statement but maybe there is more elegant solution? Yes, work with smaller buffers... Thanks again for all yours suggestions guys. Will try to implement this. Bye -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Possible bug. Clipping raster with gdalwarp.
Hi, when I clip singleband raster in GeoTiff format using shapefile as mask layer with gdalwarp I get two bands raster in output. First band is clipped raster and second - rasterized mask. Is this bug or feature? Is it really necessary to save rasterized mask? I use next command gdalwarp -q -cutline mask.shp -dstalpha -of GTiff clearcuts.tif clearcuts_clip.tif -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Merge rasters by bands
Hi all, is it possible to merge several multiband rasters in single multiband raster by bands? For example, merge two 6-bands rasters into single 12-bands raster. Using gdal_merge.py I get only 2 bands with or without -separate option. Any hints how to do this? Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Problems building GDAL 1.8.0 from sources
Hi all, I try to build GDAL 1.8.0 under Salckware 13.1 with GRASS support. As first step I build GDAL without GRASS support with this options ./configure --prefix=/usr \ --exec-prefix=$PKG/usr \ --with-pymoddir=/usr/lib/python/site-packages \ --with-libtiff=/usr/lib \ --with-geotiff=/usr/lib \ --with-geos=yes \ --with-jasper=yes \ --with-sqlite3=/usr \ --with-poppler=yes \ --with-static-proj4 \ --with-python \ --program-prefix= \ --program-suffix= \ --build=$CHOST-slackware-linux With this options GDAL uses system (external) libs and libtiff without BigTIFF LIBTIFF support: external (BigTIFF=no) But ./configure script can't find SQLite checking for sqlite3_open in -lsqlite3... no checking for SQLite3 library = 3.0.0... disabled SQLite 3.6.23.1 installed and all libs and includes available. I think this is because libdl - the dynamic linker library don't automatically supported on Slackware. This is first issue. Ok, I build GDAL without SQLite, then I build and install GRASS 6.4.1RC1. Now when I try to configure GDAL with GRASS as ./configure --prefix=/usr \ --exec-prefix=$PKG/usr \ --with-pymoddir=/usr/lib/python/site-packages \ --with-libtiff=/usr/lib \ --with-geotiff=/usr/lib \ --with-geos=yes \ --with-jasper=yes \ --with-sqlite3=/usr \ --with-poppler=yes \ --with-grass=/usr/grass-6.4.1RC1 \ --with-static-proj4 \ --with-python \ --program-prefix= \ --program-suffix= \ --build=$CHOST-slackware-linux I get strange output about libtiff support LIBTIFF support: external (BigTIFF=yes) Strange, previous build have BigTIFF=no... Also building with this config failed with errors geotiff.cpp: In member function ‘virtual CPLErr GTiffRasterBand::SetColorTable(GDALColorTable*)’: geotiff.cpp:1383: error: ‘TIFFUnsetField’ was not declared in this scope geotiff.cpp: In member function ‘void GTiffDataset::WriteGeoTIFFInfo()’: geotiff.cpp:4183: error: ‘TIFFUnsetField’ was not declared in this scope geotiff.cpp: In static member function ‘static int GTiffDataset::WriteMetadata(GDALDataset*, TIFF*, int, const char*, const char*, char**, int)’: geotiff.cpp:4635: error: ‘TIFFUnsetField’ was not declared in this scope make[2]: *** [../o/geotiff.lo] Error 1 make[2]: Leaving directory `/tmp/txz/gdal-1.8.0/frmts/gtiff' make[1]: *** [gtiff-install-obj] Error 2 make[1]: Leaving directory `/tmp/txz/gdal-1.8.0/frmts' make: *** [frmts-target] Error 2 Any thoughts why this happens and how can be fixed? 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
Thanks Mateusz, just tested and all works fine -- Alexander Bruy ___ 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] RFC 29 - OGR Ignored Fields Applied in Trunk
Thanks for this Frank We used this functionality in our project and were forced to compile GDAL with this patch. Also we add this capability to the SQLite driver. I've upload our patch for it this evening 2010/10/19 Frank Warmerdam warmer...@pobox.com: Folks, With minor adjustments and the addition of swig bindings and a test suite I have committed Martin's changes for RFC 29 - OGR field ignoring. Testing is welcome. We could of course benefit from implementation of this capability in many other drivers - especially RDBMS drivers, and it would be nice to have some applications taking advantage of it. http://trac.osgeo.org/gdal/ticket/3716 http://trac.osgeo.org/gdal/wiki/rfc29_desired_fields 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 -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Accessing in-memory shapefiles and db
Hi all, I've develop a small GIS tool and use GDAL/OGR to work with raster and vector data. Is it possible to access shapefiles which loaded in memory with OGR? And same things with SpatiaLite databases? I know that SQLIte/SpatiaLite can create in-memory DB. But is it possible to access such databases with OGR? Small example or any hinst will be very useful Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Problem with pixel/line coordinates
Thanks for your explanation, Frank! Now it works fine On Thu, 17 Jun 2010 13:59:46 -0400 Frank Warmerdam warmer...@pobox.com wrote: Alexander, I believe the addition of 0.5 in the above is incorrect. In the simple, non-rotated case, all values from geoTransform[0] to geoTransform[0] + geoTransform[1] would be on the 1st pixel (ie. pX = 0). As you have coded it, when you are half way across the pixel you are switching into the next one. Keep in mind that (geoTransform[0],geoTransform[3]) is the top left corner of the top left pixel - not the center. def pixelToMap( pX, pY, geoTransform ): mX, mY = applyGeoTransform( pX, pY, geoTransform ) return mX, mY Conversely, here if pX and pY are coming in as integer rather than floating point values, then you will likely want to add half a pixel before transforming so that you get the geoeferenced location of the center of the pixel rather than the upper left corner. -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Problem with pixel/line coordinates
Hi, I use next code to translate real coordinates (lat/lon) to the pixel/line coordinates def mapToPixel( mX, mY, geoTransform ): if geoTransform[ 2 ] + geoTransform[ 4 ] == 0: pX = ( mX - geoTransform[ 0 ] ) / geoTransform[ 1 ] pY = ( mY - geoTransform[ 3 ] ) / geoTransform[ 5 ] else: pX, pY = applyGeoTransform( mX, mY, invertGeoTransform( geoTransform ) ) return int( pX + 0.5 ), int( pY + 0.5 ) def pixelToMap( pX, pY, geoTransform ): mX, mY = applyGeoTransform( pX, pY, geoTransform ) return mX, mY def applyGeoTransform( inX, inY, geoTransform ): outX = geoTransform[ 0 ] + inX * geoTransform[ 1 ] + inY * geoTransform[ 2 ] outY = geoTransform[ 3 ] + inX * geoTransform[ 4 ] + inY * geoTransform[ 5 ] return outX, outY def invertGeoTransform( geoTransform ): det = geoTransform[ 1 ] * geoTransform[ 5 ] - geoTransform[ 2 ] * geoTransform[ 4 ] if abs( det ) 0.001: return invDet = 1.0 / det # compute adjoint and divide by determinate outGeoTransform = [ 0, 0, 0, 0, 0, 0 ] outGeoTransform[ 1 ] = geoTransform[ 5 ] * invDet outGeoTransform[ 4 ] = -geoTransform[ 4 ] * invDet outGeoTransform[ 2 ] = -geoTransform[ 2 ] * invDet outGeoTransfrom[ 5 ] = geoTransform[ 1 ] * invDet outGeoTransform[ 0 ] = ( geoTransform[ 2 ] * geoTransform[ 3 ] - geoTransform[ 0 ] * geoTransform[ 5 ] ) * invDet outGeoTransform[ 3 ] = ( -geoTransform[ 1 ] * geoTransform[ 3 ] + geoTransform[ 0 ] * geoTransform[ 4 ] ) * invDet return outGeoTransform When I test it on large raster (in AIG/Arc/Info Binary Grid format) I get next problem: raster value in some points reported by script is differ from value reported by Info tool in QGIS and/or ArcGIS. But for another points reported raster value is correct. As I understand, when point is near pixel border, MapToPixel function return wrong result - row and column of neighbour pixel. For example, if point is near right side of the pixel I get column and row for right neighbour pixel not for pixel I need. Same problem when point is near top side of pixel, I get colunm and row for top neighbour pixel. Is it possible to get correct result for all points? May be there is a mistake in code. Any help and suggestions are welcome. I can upload sample data if necessary but they are large (~330 Mb) I'm use GDAL 1.6.3 with Python 2.5.2 under Linux. Thanks! -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Handling CPG (encoding) file
Hi, GIS-Lab community has worked on this and here is our results. We develop a patch which add two methods GetEncoding and SetEncoding. At this time only GetEncoding for shapefiles is realised. With this methods programmer can get encoding programmatically. Patch uploaded to the bugtracker: http://trac.osgeo.org/gdal/ticket/882 How this works: first we try to read DBF LDID and get encoding from here (even CPG file is available), if LDID is 0 or 87 (system encoding) driver will try to read encoding from CPG file. Encoding in CPG file can be written according to the DBF LanguageID standard or according to the codepages ANSI, OEM etc. GetEncoding method will handle all this variants automatically. We'd appreciate if someone could have a look and sign off or comment on the proposed patch. On Sun, 23 May 2010 19:11:46 +0300 Alexander Bruy alexander.b...@gmail.com wrote: Hi devs, few days ago in qgis-developer mailing list we discuss codepage handling in QGIS and GDAL http://osgeo-org.1803224.n2.nabble.com/adding-encoding-file-td5075541.html. Also there is post in this list with related problem http://osgeo-org.1803224.n2.nabble.com/gdal-dev-DBF-unicode-problem-td5075730.html Seems than shapelib supports codepages in CPG file and LDID both for reading and creating DBFs. But OGR doesn't seem to use or expose that functionality. Our community want to add to the OGR shapefile driver ability to read/write CPG file and ready to work on this and create a patch. Any hints and suggestions are welcome. We really need this functionality Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Handling CPG (encoding) file
Hi devs, few days ago in qgis-developer mailing list we discuss codepage handling in QGIS and GDAL http://osgeo-org.1803224.n2.nabble.com/adding-encoding-file-td5075541.html. Also there is post in this list with related problem http://osgeo-org.1803224.n2.nabble.com/gdal-dev-DBF-unicode-problem-td5075730.html Seems than shapelib supports codepages in CPG file and LDID both for reading and creating DBFs. But OGR doesn't seem to use or expose that functionality. Our community want to add to the OGR shapefile driver ability to read/write CPG file and ready to work on this and create a patch. Any hints and suggestions are welcome. We really need this functionality Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Create new fields in existing shapefile with Python
Hi all, when writing a script in Python for processing shapefiles I found that I can add to the shapefile new fields with the same name many times without any error or warning. But AFAIK all fields in shapefile must have unique field names. Is this a bug and I need to post a ticket? Or maybe I do something wrong? Here is code #!/usr/bin/env python # -*- coding: utf-8 -*- try: from osgeo import ogr except ImportError: import ogr import sys def createField( inLayer, fieldName ): fieldDef = ogr.FieldDefn( fieldName, ogr.OFTReal ) fieldDef.SetWidth( 12 ) fieldDef.SetPrecision( 4 ) if inLayer.CreateField( fieldDef ) != 0: print Can't create field %s % fieldDef.GetNameRef() return False return True if __name__ == '__main__': args = sys.argv[ 1: ] filePath = args[ 0 ] inShape = ogr.Open( filePath, True ) if inShape is None: print 'Unable to open shapefile', filePath sys.exit( 1 ) inLayer = inShape.GetLayer( 0 ) for i in range( 3 ): name = 'new_field' if not createField( inLayer, name ): print 'ERROR' sys.exit( 1 ) Script get from commandline path to any shapefile and adds to it 4 new fields with same name. My software GDAL/OGR 1.6.3 with Python 2.5.2 under Linux. Thanks -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Real and raster coordinates
Thanks! That what I need 2010/3/25 Pinner, Luke luke.pin...@environment.gov.au: The following module might help, use the mapToPixel() function. #---Imports try: from osgeo import gdal except ImportError: import gdal def mapToPixel(mx,my,gt): ''' Convert map to pixel coordinates �...@param mx Input map x coordinate (double) �...@param my Input map y coordinate (double) �...@param gt Input geotransform (six doubles) �...@return px,py Output coordinates (two doubles) ''' if gt[2]+gt[4]==0: #Simple calc, no inversion required px = (mx - gt[0]) / gt[1] py = (my - gt[3]) / gt[5] else: px,py=ApplyGeoTransform(mx,my,InvGeoTransform(gt)) return int(px+0.5),int(py+0.5) def pixelToMap(px,py,gt): ''' Convert pixel to map coordinates �...@param px Input pixel x coordinate (double) �...@param py Input pixel y coordinate (double) �...@param gt Input geotransform (six doubles) �...@return mx,my Output coordinates (two doubles) ''' mx,my=ApplyGeoTransform(px,py,gt) return mx,my def ApplyGeoTransform(inx,iny,gt): ''' Apply a geotransform �...@param inx Input x coordinate (double) �...@param iny Input y coordinate (double) �...@param gt Input geotransform (six doubles) �...@return outx,outy Output coordinates (two doubles) ''' outx = gt[0] + inx*gt[1] + iny*gt[2] outy = gt[3] + inx*gt[4] + iny*gt[5] return (outx,outy) def InvGeoTransform(gt_in): -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Real and raster coordinates
Hi, I have a georeferenced raster and coordinates of some point (lat, lon) inside raster. Is it possible to convert this lat/lon coordinates to raster coordinates (column and row)? I use GDAL 1.6.3 with Python 2.5 (installed with OSGeo4W installer). -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Raster algebra with Gdal and Python
2010/2/1 Greg Coats gregco...@mac.com: Yes, GDAL can successfully process images with 15,000 columns by 15,000 rows. I have used GDAL on a computer with only 3 GB RAM to process an image with Image Width: 260,000 Image Length: 195,000. What version of GAL are you using? Are you using GDAL 1.7.0? What operating system are you using? How much RAM is available to GDAL? Greg I'm using GDAL 1.6.3 under Linux Slackware 12.2, and I plan to upgrade to 1.7.0 today. My computer has 2 Gb RAM and 2 Gb swap, but now I can't test this error by myself because I have no large files. This error was reported by users, they test code on several Windows XP boxes with 4 Gb RAM. I get large rasters today and try to reproduce this error. -- Alexander Bruy mailto: alexander.b...@gmail.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
