[gdal-dev] GDAL 1.9 Status
Hello, Is there any rough estimate on when GDAL/OGR 1.9 will be released ? Thanks, Silas ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Runtime error on executing OGRRegisterAll() function.
Hi, This has started happening to me also. When i call any function from gdal library i get this error: Unhandled exception at 0x7264082a (msvcr90d.dll) in test.exe: 0xC005: Access violation reading location 0x6973762f. Error happens at iosfwd file at this function static int __CLRCALL_OR_CDECL compare(_In_count_(_Count) const _Elem *_First1, _In_count_(_Count) const _Elem *_First2, size_t _Count){ // compare [_First1, _First1 + _Count) with [_First2, ...) // _DEBUG_POINTER(_First1); // _DEBUG_POINTER(_First2); return (::memcmp(_First1, _First2, _Count)); i get error at this line } I also built gdal with /MDd option but i still get this error. Any help will be appreciated. Thanks in advance.. -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/Runtime-error-on-executing-OGRRegisterAll-function-tp4453903p7096855.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] re-gridding to a coarser grid
Joaquim, nearneighbor is not really a good thing to use. Grid it with surface. What do you mean by grid it with surface? In my GDAL (1.8.1), I only have the following choices: Available resampling methods: near (default), bilinear, cubic, cubicspline, lanczos. thanks for your insight, Andreas. Joaquim Sent from my iPadiola On 14/12/2011, at 18:32, Andreas H. li...@hilboll.de wrote: Travis, yes, thanks, I had already found that in the documentation. I'm just wondering what e.g. 'bilinear' means when I go from a fine to a coarse grid? If GDAL works in terms of nodes, then I would assume 'bilinear' means interpolation, which in turn would be a very different result from working in terms of cells and then taking e.g. the average of all old cells in the new grid. I hope you see my problem / question. Thanks for you insight, Andreas. Am 14.12.2011 16:35, schrieb Travis Kirstine: Andreas, Yes gdalwarp support various resampling methods To use different resampling methods use the -r flag followed by the method eg gdalwarp -r near .. gdalwarp -r bilinear . etc... On 14 December 2011 08:26, Andreas H. li...@hilboll.de wrote: Travis, thanks for your answer! Regarding the resampling methods: Do they all just interpolate the data? I mean, when downsampling, usually I would use mean() or something similar to fill the new (coarser) grid cells. Doas gdalwarp actually do this and I'm not able to understand the documentation, or is it different? Thanks again! Andreas. Andreas, gdalwarp can be used to resample images using the -tr flag or -ts flag. For example resample 1m image to 10m using cubic resampling and target resolution' gdalwarp -r cubic -tr 10 10 input_1m.tif output_10m.tif You may have an issue determining the output resolution as I believe the resolution will be in decimal degrees rather than arc-seconds. Regards On 13 December 2011 14:52, Andreas H. li...@hilboll.de wrote: Hi, let's say I have a GeoTIFF file with a global grid in a 30 arc-second resolution. Which would be the appropriate GDAL command to spatially down-sample this file to say 0.125°? Thanks for your insight, Andreas. ___ 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 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] Warping onto an image with only GCPs
Hi, Hróbjartur Þorsteinsson explained to me a near-hidden feature of GDAL, namely to reproject one image onto the coverage of another. Given two images: imageA.tiff imageB.tiff Then an empty image can be created from image A with e.g. $ gdal_translate -ot Float32 -scale 0 0 999 999 -a_nodata 999 imageA.tiff domainA.tiff And then imageB can be warped onto the extent of this image simply with: $ gdalwarp imageB.tiff domainA.tiff domainA.tiff will then contain the data from imageB, exactly co-located with imageA. The problem is that this only works when imageA contains a geotransform, and not when it contains only GCPs. So my question is: Is there any workaround to warp data from one image onto another image which is geolocated by GCPs only? With the -to option of gdalwarp it is possible to pass some parameters to GDALCreateGenImgProjTransformer2(). In the documentation of this function I find the following option: METHOD: may have a value which is one of GEOTRANSFORM, GCP_POLYNOMIAL, GCP_TPS, GEOLOC_ARRAY, RPC to force only one geolocation method to be considered on the source dataset. http://www.gdal.org/gdal__alg_8h.html#94cd172f78dbc41d6f407d662914f2e3 But there is no corresponding option to force use of GCP_POLYNOMIAL for the *destination* dataset. I guess this is a bad sign? Thank you for any help! Knut-Frode ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Runtime error on executing OGRRegisterAll() function.
On 15 December 2011 11:49, mrym meryem...@gmail.com wrote: Hi, This has started happening to me also. When i call any function from gdal library i get this error: Unhandled exception at 0x7264082a (msvcr90d.dll) in test.exe: 0xC005: Access violation reading location 0x6973762f. Error happens at iosfwd file at this function static int __CLRCALL_OR_CDECL compare(_In_count_(_Count) const _Elem *_First1, _In_count_(_Count) const _Elem *_First2, size_t _Count){ // compare [_First1, _First1 + _Count) with [_First2, ...) // _DEBUG_POINTER(_First1); // _DEBUG_POINTER(_First2); return (::memcmp(_First1, _First2, _Count)); i get error at this line } I'm suspicious it is crashing near std::map and the crash is due to mixed variants of C/C++ runt-ime indeed. It would be good if you can post a minimal complete test program crashing with this error and don't forget to add compilation flags you used to build it. #include .your headers here.. int main() { your code here } I also built gdal with /MDd option but i still get this error. Any help will be appreciated. Thanks in advance.. Did you just changed flags and build or was it clean and complete *rebuild*? Best regards, -- Mateusz Loskot, http://mateusz.loskot.net Charter Member of OSGeo, http://osgeo.org Member of ACCU, http://accu.org ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Runtime error on executing OGRRegisterAll() function.
Here is my command line: /Od /I E:\BELGIS\gdal-1.7.2\libecwj2-3.3\Source\include /I C:\mrym\warmerda\bld\include /D WIN32 /D _DEBUG /D _WINDOWS /D _USRDLL /D BELGDAL_EXPORTS /D _WINDLL /D _UNICODE /D UNICODE /Gm /EHsc /RTC1 /MDd /FoDebug\\ /FdDebug\vc90.pdb /W3 /nologo /c /ZI /Gd /TP /errorReport:prompt *Here is my test code:* #include iostream #include windows.h #include tchar.h typedef void (*ProjectPointCall)(); int main(){ set_PGdal_ProjType_In _set_PGdal_ProjType_In; set_PGdal_ProjType_Out _set_PGdal_ProjType_Out; ProjectPointCall _ProjectPointCall; HINSTANCE hInstLibrary = LoadLibrary(_T(belGdal.dll)); if (hInstLibrary){ _ProjectPointCall=(ProjectPointCall)GetProcAddress(hInstLibrary, ProjectPointCall); _ProjectPointCall(); FreeLibrary(hInstLibrary); } else{ std::cout DLL Failed To Load! std::endl; } return 0; } *This is _ProjectPointCall():* void __cdecl ProjectPointCall(){ double *aXOut= new double; double *aYOut= new double; Gdal_ProjectPoint( 0.0, 0.0, p_ProjType_In, p_ProjType_Out, aXOut, aYOut ); } *p_ProjType_In, p_ProjType_Out are static variables. This is Gdal_ProjectPoint function:* bool __stdcall Gdal_ProjectPoint(double aXIn, double aYIn, const PGdal_ProjType aProjIn, const PGdal_ProjType aProjOut, double* aXOut, double* aYOut ){ bool changeProjSrc, changeProjDest; int datumCount = g_EPSGDatums.GetSize(); GDALAllRegister(); if(datumCount == 0) SetupEPSGDatums(); } When i debug it i get error at GDALAllRegister(). I've tried debugging other functions in my dll and all gdal function calls gave me the same error. And i have changed flags and built it in this way nmake -f makefile.vc clean nmake -f makefile.vc MSVC_VER=1500 DEBUG=1 nmake /f makefile.vc install nmake /f makefile.vc devinstall Thanks a lot.. -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/Runtime-error-on-executing-OGRRegisterAll-function-tp4453903p7097491.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] Slope in gdaldem
Hi All, A colleague tipped me off to an article in Forestry Source entitled When GIS Slope Isn't What You Think showing differences between field measured slope values and GIS results. The differences can be substantial. The short and long of it is that local hydrologic slope is closer to what's used in the field by foresters than the integrated slope from either Horn or Zevenbergen-Thorne. And so, I have a feature request. I would like to add hydrologic slope to gdaldem. This is max local slope rather than integrated 9 cell slope like the existing code: The existing code for Horn, e.g. is: dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData-ewres; dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData-nsres; I'd like something more like this (adapted from roughness): float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } return pafLocalMax - afWin[4]; } Any critques? Reasons to not add? Best, Steve Stephen Mather GIS Manager Cleveland Metroparks 4101 Fulton Pkwy Cleveland, OH 44144 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Runtime error on executing OGRRegisterAll() function.
On 15 December 2011 15:33, mrym meryem...@gmail.com wrote: And i have changed flags and built it in this way nmake -f makefile.vc clean nmake -f makefile.vc MSVC_VER=1500 DEBUG=1 Edit nmake.opt In debug version of OPTFLAGS flags, replace /MD with /MDd Rebuild GDAL Rebuild your library and program. Check if it .helps. Best regards, -- Mateusz Loskot, http://mateusz.loskot.net Charter Member of OSGeo, http://osgeo.org Member of ACCU, http://accu.org ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_translate to TIFF -- color mismatch
Hi all, I am trying to use gdal_translate to convert a .grd file (Northwood/VerticalMapper format) to TIFF. However, the colors in the translated file don't seem to be interpolated correctly based on the color scheme, for values near the end of the range. Here is a link to the original GRD, as well as the exported TIFFs (small download, only 200K): http://dl.dropbox.com/u/53500018/GdalQuestion.zip The GdalTranslateExport.tif file was produced with the command gdal_translate -b 1 -b 2 -b 3 C:\source\Ottawa_east_DEM30m_trim_copy.grd C:\temp\GdalTranslateExport.tif. Notice the top and bottom areas of the image: they are yellow. For comparison, take a look at the NativeExport.tif file in the same package. The top and bottom areas are orange, as they should be. The NativeExport.tif file was produced by Pitney Bowes MapInfo. Here is the color scheme for the GRD: http://www.smugmug.com/photos/i-bRZ879w/0/X3/i-bRZ879w-X3.jpg Notice that the thresholds are 70 for yellow, and 161 for orange. Well, the values at the top and bottom are in the high 130s/low 140s. So they should be interpolated close to orange, not yellow. Therefore the NativeExport.tif is correct, and the gdal_translate export is incorrect. I have also tried changing the thresholds, and as soon as I change the orange threshold to be below the maximum actual value in the GRD (let's say 144), the gdal_translate export shows the orange shades. So, it seems there might be a bug in gdal_translate where color values are not interpolated correctly for thresholds that fall outside of the range of actual values? Or is there a way to make gdal_translate behave the way that I want with some different command line options? Regards, Mladen Gavrilovic P.S. For convenience, here is a gallery with the TIFFs converted to JPEG (difference is still obvious), and the color scheme: http://smu.gs/vMEemF P.P.S Version of GDAL is 1.8.1 This message was sent using Distributel Webmail. Ce message a été envoyé à partir de la Messagerie Web Distributel. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] extra spaces when writing to AAIGRID format
Hi, I'm having problems with GDAL inserting leading spaces when writing from a numpy array to ESRI ASCII Grid format. The problem occurs whenever there are integers in the array (even though the array itself is floating point). E.g., if I try to write a 3x3 grid that ArcGIS can read, it needs to look something like this: ncols3 nrows3 xllcorner100. yllcorner100. cellsize 0.1000 1.0 1.4809201955795288 1.4809201955795288 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10.00 Instead, the ascii file written by GDAL looks like this (which can't be read by ArcGIS because of extra leading spaces in the first line of the data): ncols5 nrows5 xllcorner100. yllcorner100. cellsize 0.1000 1 1.4809201955795288 1.4809201955795288 0.25947147607803345 0.71342569589614868 2.810183288701 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10 0.74458205699920654 Here is the code I have been trying to use: format = MEM driver = gdal.GetDriverByName( format ) dst_ds = driver.Create( file, len(data[0]), len(data),1,gdal.GDT_Float32) dst_ds.SetGeoTransform([xllcorner, cellsize, 0, ull, 0, cellsize]) dst_ds.GetRasterBand(1).WriteArray(data) format = 'AAIGrid' driver = gdal.GetDriverByName(format) dst_ds_new = driver.CreateCopy(file, dst_ds) dst_ds = None I've seen this problem on various versions of python going back at least 3 years, so have relied on the much slower numpy.savetxt command instead of using GDAL. Any thoughts on how to make this work? Thanks! Brad -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/extra-spaces-when-writing-to-AAIGRID-format-tp7097852p7097852.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] re-gridding to a coarser grid
Andreas I think thats referring to my suggestion of GMT to do part of the work On Thu, 2011-12-15 at 14:00 +0100, Andreas H. wrote: Joaquim, nearneighbor is not really a good thing to use. Grid it with surface. What do you mean by grid it with surface? In my GDAL (1.8.1), I only have the following choices: Available resampling methods: near (default), bilinear, cubic, cubicspline, lanczos. thanks for your insight, Andreas. Joaquim Sent from my iPadiola On 14/12/2011, at 18:32, Andreas H. li...@hilboll.de wrote: Travis, yes, thanks, I had already found that in the documentation. I'm just wondering what e.g. 'bilinear' means when I go from a fine to a coarse grid? If GDAL works in terms of nodes, then I would assume 'bilinear' means interpolation, which in turn would be a very different result from working in terms of cells and then taking e.g. the average of all old cells in the new grid. I hope you see my problem / question. Thanks for you insight, Andreas. Am 14.12.2011 16:35, schrieb Travis Kirstine: Andreas, Yes gdalwarp support various resampling methods To use different resampling methods use the -r flag followed by the method eg gdalwarp -r near .. gdalwarp -r bilinear . etc... On 14 December 2011 08:26, Andreas H. li...@hilboll.de wrote: Travis, thanks for your answer! Regarding the resampling methods: Do they all just interpolate the data? I mean, when downsampling, usually I would use mean() or something similar to fill the new (coarser) grid cells. Doas gdalwarp actually do this and I'm not able to understand the documentation, or is it different? Thanks again! Andreas. Andreas, gdalwarp can be used to resample images using the -tr flag or -ts flag. For example resample 1m image to 10m using cubic resampling and target resolution' gdalwarp -r cubic -tr 10 10 input_1m.tif output_10m.tif You may have an issue determining the output resolution as I believe the resolution will be in decimal degrees rather than arc-seconds. Regards On 13 December 2011 14:52, Andreas H. li...@hilboll.de wrote: Hi, let's say I have a GeoTIFF file with a global grid in a 30 arc-second resolution. Which would be the appropriate GDAL command to spatially down-sample this file to say 0.125°? Thanks for your insight, Andreas. ___ 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 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 mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] re-gridding to a coarser grid
Sorry for the noise that I introduced by mixing the GDAL and GMT's lists. I already replied to Andreas only because I though my message went only to him, but since it didn't here is a copy of it. Joaquim Andreas, I'm very sorry but I sent that message to you by mistake as it was meant as a reply to another question in the GMT list. But since we are at it, your case is also probably better suite for GMT than GDAL. In GMT you have 'grdsample' to resample your grid or, if you really want to do things the correct way to avoid aliasing since you are going from high to lower resolution, 'grdfilter'. Note that with GMT5 built with GDAL support you can read you grid directly but you are more limited on output formats (mostly netCDF). Andreas I think thats referring to my suggestion of GMT to do part of the work On Thu, 2011-12-15 at 14:00 +0100, Andreas H. wrote: Joaquim, nearneighbor is not really a good thing to use. Grid it with surface. What do you mean by grid it with surface? In my GDAL (1.8.1), I only have the following choices: Available resampling methods: near (default), bilinear, cubic, cubicspline, lanczos. thanks for your insight, Andreas. Joaquim Sent from my iPadiola On 14/12/2011, at 18:32, Andreas H.li...@hilboll.de wrote: Travis, yes, thanks, I had already found that in the documentation. I'm just wondering what e.g. 'bilinear' means when I go from a fine to a coarse grid? If GDAL works in terms of nodes, then I would assume 'bilinear' means interpolation, which in turn would be a very different result from working in terms of cells and then taking e.g. the average of all old cells in the new grid. I hope you see my problem / question. Thanks for you insight, Andreas. Am 14.12.2011 16:35, schrieb Travis Kirstine: Andreas, Yes gdalwarp support various resampling methods To use different resampling methods use the -r flag followed by the method eg gdalwarp -r near .. gdalwarp -r bilinear . etc... On 14 December 2011 08:26, Andreas H.li...@hilboll.de wrote: Travis, thanks for your answer! Regarding the resampling methods: Do they all just interpolate the data? I mean, when downsampling, usually I would use mean() or something similar to fill the new (coarser) grid cells. Doas gdalwarp actually do this and I'm not able to understand the documentation, or is it different? Thanks again! Andreas. Andreas, gdalwarp can be used to resample images using the -tr flag or -ts flag. For example resample 1m image to 10m using cubic resampling and target resolution' gdalwarp -r cubic -tr 10 10 input_1m.tif output_10m.tif You may have an issue determining the output resolution as I believe the resolution will be in decimal degrees rather than arc-seconds. Regards On 13 December 2011 14:52, Andreas H.li...@hilboll.de wrote: Hi, let's say I have a GeoTIFF file with a global grid in a 30 arc-second resolution. Which would be the appropriate GDAL command to spatially down-sample this file to say 0.125°? Thanks for your insight, Andreas. ___ 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 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 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
Re: [gdal-dev] re-gridding to a coarser grid
If you need complete control, and your data are not too massive in number, I have had good results pushing pixel coordinates as points into postgis and then using spatial queries to aggregate in various ways (eg averages a hexagonal grid). It seems round-about, but it works. I've actually done it on some pretty massive grids. THK On Wed, Dec 14, 2011 at 7:26 AM, Andreas H. li...@hilboll.de wrote: Travis, thanks for your answer! Regarding the resampling methods: Do they all just interpolate the data? I mean, when downsampling, usually I would use mean() or something similar to fill the new (coarser) grid cells. Doas gdalwarp actually do this and I'm not able to understand the documentation, or is it different? Thanks again! Andreas. Andreas, gdalwarp can be used to resample images using the -tr flag or -ts flag. For example resample 1m image to 10m using cubic resampling and target resolution' gdalwarp -r cubic -tr 10 10 input_1m.tif output_10m.tif You may have an issue determining the output resolution as I believe the resolution will be in decimal degrees rather than arc-seconds. Regards On 13 December 2011 14:52, Andreas H. li...@hilboll.de wrote: Hi, let's say I have a GeoTIFF file with a global grid in a 30 arc-second resolution. Which would be the appropriate GDAL command to spatially down-sample this file to say 0.125°? Thanks for your insight, Andreas. ___ 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 -- Timothy H. Keitt http://www.keittlab.org/ ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Slope in gdaldem
I already see a flaw-- the function should be more like: float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells const double radiansToDegrees = 180.0 / M_PI; GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } return atan(sqrt(pafLocalMax - afWin[4]) / (2*psData-scale)) * radiansToDegrees; } On Thu, Dec 15, 2011 at 11:49 AM, Stephen Mather mather.step...@gmail.comwrote: Hi All, A colleague tipped me off to an article in Forestry Source entitled When GIS Slope Isn't What You Think showing differences between field measured slope values and GIS results. The differences can be substantial. The short and long of it is that local hydrologic slope is closer to what's used in the field by foresters than the integrated slope from either Horn or Zevenbergen-Thorne. And so, I have a feature request. I would like to add hydrologic slope to gdaldem. This is max local slope rather than integrated 9 cell slope like the existing code: The existing code for Horn, e.g. is: dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData-ewres; dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData-nsres; I'd like something more like this (adapted from roughness): float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } return pafLocalMax - afWin[4]; } Any critques? Reasons to not add? Best, Steve Stephen Mather GIS Manager Cleveland Metroparks 4101 Fulton Pkwy Cleveland, OH 44144 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Slope in gdaldem
Apologies for the flurry, one more change: float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells const double radiansToDegrees = 180.0 / M_PI; GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData; double key; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } key = pafLocalMax - afWin[4]; if (psData-slopeFormat == 1) return (float) (atan(sqrt(key) / (2*psData-scale)) * radiansToDegrees); else return (float) (100*(sqrt(key) / (2*psData-scale))); } Best, Steve ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Slope in gdaldem
Stephen, It'd be great to see local hydrologic slope in gdaldem. I've seen a similar discrepancy between field slope and most GIS-derived slopes. The local max hydrologic slope seems to correspond much better to our perceived slope on the ground. The integrated slope used in most GIS, however, better matches the overall planar slope over a large extent. So they both have their place depending on scale; local hydrologic max slope for the micro, integrated planar slope for the macro. However, I don't think the algorithm that you describe below will work. It's currently looking for the *maximum* local elevation and computing slope to the center cell. Instead it should calculate the slope in the direction of hydrologic flow; i.e. the slope from the center cell to surrounding cell with the *minimum* elevation. - matt On Thu, Dec 15, 2011 at 10:29 AM, Stephen Mather mather.step...@gmail.com wrote: I already see a flaw-- the function should be more like: float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells const double radiansToDegrees = 180.0 / M_PI; GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } return atan(sqrt(pafLocalMax - afWin[4]) / (2*psData-scale)) * radiansToDegrees; } On Thu, Dec 15, 2011 at 11:49 AM, Stephen Mather mather.step...@gmail.com wrote: Hi All, A colleague tipped me off to an article in Forestry Source entitled When GIS Slope Isn't What You Think showing differences between field measured slope values and GIS results. The differences can be substantial. The short and long of it is that local hydrologic slope is closer to what's used in the field by foresters than the integrated slope from either Horn or Zevenbergen-Thorne. And so, I have a feature request. I would like to add hydrologic slope to gdaldem. This is max local slope rather than integrated 9 cell slope like the existing code: The existing code for Horn, e.g. is: dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData-ewres; dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData-nsres; I'd like something more like this (adapted from roughness): float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } return pafLocalMax - afWin[4]; } Any critques? Reasons to not add? Best, Steve Stephen Mather GIS Manager Cleveland Metroparks 4101 Fulton Pkwy Cleveland, OH 44144 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- ·´¯`·.¸. , . .·´¯`·.. º`·.¸¸.·´¯`·.¸.·´¯`·...¸º The best way to predict the future is to invent it. -- Alan Kay Matthew T. Perry http://www.perrygeo.net http://viedevelo.wordpress.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] re-gridding to a coarser grid
For what it's worth, What do you consider large and what kind of method for downscaling you wish to use? I do this in GRASS using the r.neighbors command, http://grass.fbk.eu/grass64/manuals/html64_user/r.neighbors.html. I've done an averaging method to convert a 20 ft elevation grid ( statewide coverage) of approximately 7 billion cells to an averaged 60 ft grid ( 777 million cells) by performing a 3x3 neighbor hood analysis on the 20 ft grid, setting the GRASS workspace resolution to 60 ft and multiplying the 20ft grid created by the 3x3 neighborhood analysis by 1 to create a new grid ( http://grass.fbk.eu/grass64/manuals/html64_user/r.mapcalculator.html) . GRASS will take the value of the center averaged 20 ft cell in the 60ft block as the value of the cell in the new 60 ft raster. options for r.neigbors calculations are : average,median,mode,minimum,maximum,range,stddev,sum,variance,diversity,interspersion Doug Doug Newcomb USFWS Raleigh, NC 919-856-4520 ext. 14 doug_newc...@fws.gov - The opinions I express are my own and are not representative of the official policy of the U.S.Fish and Wildlife Service or Dept. of the Interior. Life is too short for undocumented, proprietary data formats. Tim Keitt tke...@gmail.com Sent by: gdal-dev-boun...@lists.osgeo.org 12/15/2011 01:17 PM To Andreas H. li...@hilboll.de cc gdal-dev@lists.osgeo.org Subject Re: [gdal-dev] re-gridding to a coarser grid If you need complete control, and your data are not too massive in number, I have had good results pushing pixel coordinates as points into postgis and then using spatial queries to aggregate in various ways (eg averages a hexagonal grid). It seems round-about, but it works. I've actually done it on some pretty massive grids. THK On Wed, Dec 14, 2011 at 7:26 AM, Andreas H. li...@hilboll.de wrote: Travis, thanks for your answer! Regarding the resampling methods: Do they all just interpolate the data? I mean, when downsampling, usually I would use mean() or something similar to fill the new (coarser) grid cells. Doas gdalwarp actually do this and I'm not able to understand the documentation, or is it different? Thanks again! Andreas. Andreas, gdalwarp can be used to resample images using the -tr flag or -ts flag. For example resample 1m image to 10m using cubic resampling and target resolution' gdalwarp -r cubic -tr 10 10 input_1m.tif output_10m.tif You may have an issue determining the output resolution as I believe the resolution will be in decimal degrees rather than arc-seconds. Regards On 13 December 2011 14:52, Andreas H. li...@hilboll.de wrote: Hi, let's say I have a GeoTIFF file with a global grid in a 30 arc-second resolution. Which would be the appropriate GDAL command to spatially down-sample this file to say 0.125°? Thanks for your insight, Andreas. ___ 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 -- Timothy H. Keitt http://www.keittlab.org/ ___ 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
Re: [gdal-dev] extra spaces when writing to AAIGRID format
Le jeudi 15 décembre 2011 18:21:11, bhmcrae a écrit : Hi, I'm having problems with GDAL inserting leading spaces when writing from a numpy array to ESRI ASCII Grid format. The problem occurs whenever there are integers in the array (even though the array itself is floating point). E.g., if I try to write a 3x3 grid that ArcGIS can read, it needs to look something like this: ncols3 nrows3 xllcorner100. yllcorner100. cellsize 0.1000 1.0 1.4809201955795288 1.4809201955795288 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10.00 Instead, the ascii file written by GDAL looks like this (which can't be read by ArcGIS because of extra leading spaces in the first line of the data): ncols5 nrows5 xllcorner100. yllcorner100. cellsize 0.1000 1 1.4809201955795288 1.4809201955795288 0.25947147607803345 0.71342569589614868 2.810183288701 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10 0.74458205699920654 Here is the code I have been trying to use: format = MEM driver = gdal.GetDriverByName( format ) dst_ds = driver.Create( file, len(data[0]), len(data),1,gdal.GDT_Float32) dst_ds.SetGeoTransform([xllcorner, cellsize, 0, ull, 0, cellsize]) dst_ds.GetRasterBand(1).WriteArray(data) format = 'AAIGrid' driver = gdal.GetDriverByName(format) dst_ds_new = driver.CreateCopy(file, dst_ds) dst_ds = None I've seen this problem on various versions of python going back at least 3 years, so have relied on the much slower numpy.savetxt command instead of using GDAL. Any thoughts on how to make this work? Brad, What iI can say for sure is that it is not related at all with which Python version you are using. This is specific to the AAIGRD driver. You would get the same results with gdal_translate. It would be great if you could identify more precisely the cause for the failed opening. Is it : 1) the presence of spaces at the beginning of lines (the first number of the line being an integer or a float), 2) the presence of spaces at the beginning of lines starting with an integer, 3) integer values not being written as floating point numbers For 1), 2) and 3), I've attached test1.asc, test2.asc and test3.asc that should test the hypothesis. I've also attached test0.asc for reference, which is the example that was generated by ArcGIS, and that I've altered to generate the other files. Could you open those files with ArcGIS and report which ones work ? And what is the ArcGIS version you are using ? Do other ArcGIS users that happen to read that thread have met similar issues ? Best regards, Even test0.asc Description: application/pgp-keys test1.asc Description: application/pgp-keys test2.asc Description: application/pgp-keys test3.asc Description: application/pgp-keys ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] extra spaces when writing to AAIGRID format
Brad, E.g., if I try to write a 3x3 grid that ArcGIS can read, it needs to look something like this: ncols 3 nrows 3 xllcorner 100. yllcorner 100. cellsize 0.1000 1.0 1.4809201955795288 1.4809201955795288 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10.00 Instead, the ascii file written by GDAL looks like this (which can't be read by ArcGIS because of extra leading spaces in the first line of the data): What is the ArcGIS error? I tried these files on ArcMap 9.3 and 10.0 and had no issues that I realized. (to display in a meaningful manner, I had to change the symbology which required the calculation of statistics in ArcGIS. I also changed the cellsize to .1 so I could more easily examine the pixels and values.) Below looks like a 3x4 grid that is labeled as 5x5 but even that worked for me in ArcGIS (didn't throw error but values were not where 'expected'). ncols 5 nrows 5 xllcorner 100. yllcorner 100. cellsize 0.1000 1 1.4809201955795288 1.4809201955795288 0.25947147607803345 0.71342569589614868 2.810183288701 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10 0.74458205699920654 All 4 of of Even's files worked for me too. HTH, Eli ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_translate to TIFF -- color mismatch
Le jeudi 15 décembre 2011 18:08:26, mlade...@distributel.net a écrit : Hi all, I am trying to use gdal_translate to convert a .grd file (Northwood/VerticalMapper format) to TIFF. However, the colors in the translated file don't seem to be interpolated correctly based on the color scheme, for values near the end of the range. Here is a link to the original GRD, as well as the exported TIFFs (small download, only 200K): http://dl.dropbox.com/u/53500018/GdalQuestion.zip The GdalTranslateExport.tif file was produced with the command gdal_translate -b 1 -b 2 -b 3 C:\source\Ottawa_east_DEM30m_trim_copy.grd C:\temp\GdalTranslateExport.tif. Notice the top and bottom areas of the image: they are yellow. For comparison, take a look at the NativeExport.tif file in the same package. The top and bottom areas are orange, as they should be. The NativeExport.tif file was produced by Pitney Bowes MapInfo. Here is the color scheme for the GRD: http://www.smugmug.com/photos/i-bRZ879w/0/X3/i-bRZ879w-X3.jpg Notice that the thresholds are 70 for yellow, and 161 for orange. Well, the values at the top and bottom are in the high 130s/low 140s. So they should be interpolated close to orange, not yellow. Therefore the NativeExport.tif is correct, and the gdal_translate export is incorrect. I have also tried changing the thresholds, and as soon as I change the orange threshold to be below the maximum actual value in the GRD (let's say 144), the gdal_translate export shows the orange shades. Good analysis. There was indeed a bug in the way the end of the color palette is interpolated if the maximum z value is below a threshold of the color scheme. See http://trac.osgeo.org/gdal/ticket/4395 for the fix (a simple typo) So, it seems there might be a bug in gdal_translate where color values are not interpolated correctly for thresholds that fall outside of the range of actual values? Or is there a way to make gdal_translate behave the way that I want with some different command line options? Regards, Mladen Gavrilovic P.S. For convenience, here is a gallery with the TIFFs converted to JPEG (difference is still obvious), and the color scheme: http://smu.gs/vMEemF P.P.S Version of GDAL is 1.8.1 This message was sent using Distributel Webmail. Ce message a été envoyé à partir de la Messagerie Web Distributel. ___ 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
Re: [gdal-dev] extra spaces when writing to AAIGRID format
Hi Even and Eli, Thanks much for your responses. As soon as I sent that email I realized I might eat my words on the inability to open the ASCII grids in ArcMap. This problem did crop up when we released a version of our code that used GDAL a few years back. A user reported intermittent errors reading grids that our code exported. Since then, I have assumed that we must fix the writing behavior of our GDAL code. But, I was able to open the files you sent in ArcMap 10 (as well as the grid example that I sent out). So it is not always a problem. Moreover, it may not be a problem at all with the latest versions of ArcGIS- the old problem may have even been on ArcView 3.x. We'll do some further testing and decide whether it seems risky at all. Thanks again, Brad On Thu, Dec 15, 2011 at 11:01 AM, Even Rouault even.roua...@mines-paris.org wrote: Le jeudi 15 décembre 2011 18:21:11, bhmcrae a écrit : Hi, I'm having problems with GDAL inserting leading spaces when writing from a numpy array to ESRI ASCII Grid format. The problem occurs whenever there are integers in the array (even though the array itself is floating point). E.g., if I try to write a 3x3 grid that ArcGIS can read, it needs to look something like this: ncols 3 nrows 3 xllcorner 100. yllcorner 100. cellsize 0.1000 1.0 1.4809201955795288 1.4809201955795288 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10.00 Instead, the ascii file written by GDAL looks like this (which can't be read by ArcGIS because of extra leading spaces in the first line of the data): ncols 5 nrows 5 xllcorner 100. yllcorner 100. cellsize 0.1000 1 1.4809201955795288 1.4809201955795288 0.25947147607803345 0.71342569589614868 2.810183288701 0.11773009598255157 0.15476560592651367 0.86710774898529053 0.17429724335670471 0.76037305593490601 10 0.74458205699920654 Here is the code I have been trying to use: format = MEM driver = gdal.GetDriverByName( format ) dst_ds = driver.Create( file, len(data[0]), len(data),1,gdal.GDT_Float32) dst_ds.SetGeoTransform([xllcorner, cellsize, 0, ull, 0, cellsize]) dst_ds.GetRasterBand(1).WriteArray(data) format = 'AAIGrid' driver = gdal.GetDriverByName(format) dst_ds_new = driver.CreateCopy(file, dst_ds) dst_ds = None I've seen this problem on various versions of python going back at least 3 years, so have relied on the much slower numpy.savetxt command instead of using GDAL. Any thoughts on how to make this work? Brad, What iI can say for sure is that it is not related at all with which Python version you are using. This is specific to the AAIGRD driver. You would get the same results with gdal_translate. It would be great if you could identify more precisely the cause for the failed opening. Is it : 1) the presence of spaces at the beginning of lines (the first number of the line being an integer or a float), 2) the presence of spaces at the beginning of lines starting with an integer, 3) integer values not being written as floating point numbers For 1), 2) and 3), I've attached test1.asc, test2.asc and test3.asc that should test the hypothesis. I've also attached test0.asc for reference, which is the example that was generated by ArcGIS, and that I've altered to generate the other files. Could you open those files with ArcGIS and report which ones work ? And what is the ArcGIS version you are using ? Do other ArcGIS users that happen to read that thread have met similar issues ? Best regards, Even -- -- Brad McRae The Nature Conservancy 1917 1st Avenue Seattle, WA 98101 Cell: 541-223-1170 Office: 206-436-6206 Fax: 206-343-5608 email: bmc...@tnc.org http://www.nceas.ucsb.edu/~mcrae/ ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Slope in gdaldem
I agree both have their place. Sometimes we get a more consistent or more interesting answer from our analyses than from our field measurements. As to the approach, I agree with your assessment, I was thinking backwards-- one caveat-- how best to handle local minima? Insert a null? I'm afraid I don't know how this is normally handled. Steve On Thu, Dec 15, 2011 at 1:45 PM, Matt Perry perry...@gmail.com wrote: Stephen, It'd be great to see local hydrologic slope in gdaldem. I've seen a similar discrepancy between field slope and most GIS-derived slopes. The local max hydrologic slope seems to correspond much better to our perceived slope on the ground. The integrated slope used in most GIS, however, better matches the overall planar slope over a large extent. So they both have their place depending on scale; local hydrologic max slope for the micro, integrated planar slope for the macro. However, I don't think the algorithm that you describe below will work. It's currently looking for the *maximum* local elevation and computing slope to the center cell. Instead it should calculate the slope in the direction of hydrologic flow; i.e. the slope from the center cell to surrounding cell with the *minimum* elevation. - matt On Thu, Dec 15, 2011 at 10:29 AM, Stephen Mather mather.step...@gmail.com wrote: I already see a flaw-- the function should be more like: float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells const double radiansToDegrees = 180.0 / M_PI; GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } return atan(sqrt(pafLocalMax - afWin[4]) / (2*psData-scale)) * radiansToDegrees; } On Thu, Dec 15, 2011 at 11:49 AM, Stephen Mather mather.step...@gmail.com wrote: Hi All, A colleague tipped me off to an article in Forestry Source entitled When GIS Slope Isn't What You Think showing differences between field measured slope values and GIS results. The differences can be substantial. The short and long of it is that local hydrologic slope is closer to what's used in the field by foresters than the integrated slope from either Horn or Zevenbergen-Thorne. And so, I have a feature request. I would like to add hydrologic slope to gdaldem. This is max local slope rather than integrated 9 cell slope like the existing code: The existing code for Horn, e.g. is: dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData-ewres; dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData-nsres; I'd like something more like this (adapted from roughness): float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } return pafLocalMax - afWin[4]; } Any critques? Reasons to not add? Best, Steve Stephen Mather GIS Manager Cleveland Metroparks 4101 Fulton Pkwy Cleveland, OH 44144 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- ·´¯`·.¸. , . .·´¯`·.. º`·.¸¸.·´¯`·.¸.·´¯`·...¸º The best way to predict the future is to invent it. -- Alan Kay Matthew T. Perry http://www.perrygeo.net http://viedevelo.wordpress.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_translate and TIFF resolutions(?)
Hi, I'm using gdal_translate to chip out TIFFs from a NITF using the command line: gdal_translate -ot Int32 -of Gtiff -strict -outsize 3000 3000 -srcwin 0 0 3000 3000 c:\data\image.ntf c:\data\output_image.tiff Tiffinfo (which throws an error) has a problem with the resolution field- the x and y resolutions seem to be unset. Irfanview will display the image, but also shows the resolution to be blank. When I modify the resolution within Irfanview, I'm able to load the image using any app- paint.net, photoshop, MS Paint, programmatically with OpenCV, etc. How can I force gdal_translate to set that resolution field? I'm not a TIFF expert, but it seems like it would be a standard field if tiffinfo throws an error when it's not populated. Thanks, Mike ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] re-gridding to a coarser grid
Sometimes I think GRASS can do anything... In my case I was using a 1km equal area grid over North and South America and needed to implement an interpolation scheme not pre-packaged anywhere I had looked, so I worked with large tables of points. Not highly recommended unless you have other things to do while you wait. THK On Thu, Dec 15, 2011 at 12:39 PM, doug_newc...@fws.gov wrote: For what it's worth, What do you consider large and what kind of method for downscaling you wish to use? I do this in GRASS using the r.neighbors command, http://grass.fbk.eu/grass64/manuals/html64_user/r.neighbors.html. I've done an averaging method to convert a 20 ft elevation grid ( statewide coverage) of approximately 7 billion cells to an averaged 60 ft grid ( 777 million cells) by performing a 3x3 neighbor hood analysis on the 20 ft grid, setting the GRASS workspace resolution to 60 ft and multiplying the 20ft grid created by the 3x3 neighborhood analysis by 1 to create a new grid (http://grass.fbk.eu/grass64/manuals/html64_user/r.neighbors.htmlhttp://grass.fbk.eu/grass64/manuals/html64_user/r.mapcalculator.html) . GRASS will take the value of the center averaged 20 ft cell in the 60ft block as the value of the cell in the new 60 ft raster. http://grass.fbk.eu/grass64/manuals/html64_user/r.mapcalculator.html options for r.neigbors calculations are : average,median,mode,minimum,maximum,range,stddev,sum,variance,diversity,interspersion Doug Doug Newcomb USFWS Raleigh, NC 919-856-4520 ext. 14 doug_newc...@fws.gov - The opinions I express are my own and are not representative of the official policy of the U.S.Fish and Wildlife Service or Dept. of the Interior. Life is too short for undocumented, proprietary data formats. *Tim Keitt tke...@gmail.com* Sent by: gdal-dev-boun...@lists.osgeo.org 12/15/2011 01:17 PM To Andreas H. li...@hilboll.de cc gdal-dev@lists.osgeo.org Subject Re: [gdal-dev] re-gridding to a coarser grid If you need complete control, and your data are not too massive in number, I have had good results pushing pixel coordinates as points into postgis and then using spatial queries to aggregate in various ways (eg averages a hexagonal grid). It seems round-about, but it works. I've actually done it on some pretty massive grids. THK On Wed, Dec 14, 2011 at 7:26 AM, Andreas H. li...@hilboll.de wrote: Travis, thanks for your answer! Regarding the resampling methods: Do they all just interpolate the data? I mean, when downsampling, usually I would use mean() or something similar to fill the new (coarser) grid cells. Doas gdalwarp actually do this and I'm not able to understand the documentation, or is it different? Thanks again! Andreas. Andreas, gdalwarp can be used to resample images using the -tr flag or -ts flag. For example resample 1m image to 10m using cubic resampling and target resolution' gdalwarp -r cubic -tr 10 10 input_1m.tif output_10m.tif You may have an issue determining the output resolution as I believe the resolution will be in decimal degrees rather than arc-seconds. Regards On 13 December 2011 14:52, Andreas H. li...@hilboll.de wrote: Hi, let's say I have a GeoTIFF file with a global grid in a 30 arc-second resolution. Which would be the appropriate GDAL command to spatially down-sample this file to say 0.125°? Thanks for your insight, Andreas. ___ 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 http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Timothy H. Keitt http://lists.osgeo.org/mailman/listinfo/gdal-dev http://www.keittlab.org/ ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://www.keittlab.org/ http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Timothy H. Keitt http://www.keittlab.org/ ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_translate and TIFF resolutions(?)
Mike, Can you provide the output of tiffinfo, and tiffdump on the file produced by gdal_translate? GDAL does have support for setting the resolution fields if particular source metadata is provided but normally that metadata is not available and no resolution fields should be written at all. I'd also like to know what version of gdal you are using. Best regards, Frank On Thu, Dec 15, 2011 at 12:50 PM, Mike O'Malley mike.omal...@nextcentury.com wrote: Hi, I’m using gdal_translate to chip out TIFFs from a NITF using the command line: gdal_translate -ot Int32 -of Gtiff -strict -outsize 3000 3000 –srcwin 0 0 3000 3000 c:\data\image.ntf c:\data\output_image.tiff Tiffinfo (which throws an error) has a problem with the resolution field- the x and y resolutions seem to be unset. Irfanview will display the image, but also shows the resolution to be blank. When I modify the resolution within Irfanview, I’m able to load the image using any app- paint.net, photoshop, MS Paint, programmatically with OpenCV, etc. How can I force gdal_translate to set that resolution field? I’m not a TIFF expert, but it seems like it would be a standard field if tiffinfo throws an error when it’s not populated. Thanks, Mike ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- ---+-- 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 Software Developer ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Fwd: [gdal-dev] gdal_translate and TIFF resolutions(?)
Forwarding private response back to the group. -- Forwarded message -- From: Frank Warmerdam warmer...@pobox.com Date: Thu, Dec 15, 2011 at 1:56 PM Subject: Re: [gdal-dev] gdal_translate and TIFF resolutions(?) To: Mike O'Malley mike.omal...@nextcentury.com On Thu, Dec 15, 2011 at 1:41 PM, Mike O'Malley mike.omal...@nextcentury.com wrote: Frank, Sure thing! I've attached the output from tiffinfo and tiffdump. Tiffinfo noted a few unknown fields with high numeric tags (34k-ish), but I don't know that those are relevant. I'm using GDAL 1.8.1. What source metadata does GDAL require to set the resolution fields? I don't have control over the NITFs, but if there's some way I can inject that metadata into the process I'd be happy to look into it. Or after the process, for that matter. Mike, The tiffinfo messages about unknown tags are not serious. Looking at the generated file I don't see any problems with it as far as resolution goes. It just doesn't have any resolution tags at all which is perfectly legitimate. If some apps can't consume files without resolution set then I think they are the issue. However, you can set it yourself like: gdal_translate ... \ -mo TIFFTAG_XRESOLUTION=80 \ -mo TIFFTAG_YRESOLUTION=80 \ -mo TIFFTAG_RESOLUTIONUNIT=2 #define RESUNIT_NONE 1 /* no meaningful units */ #define RESUNIT_INCH 2 /* english */ #define RESUNIT_CENTIMETER 3 /* metric */ The above would be 80 dpi. The tags are mentioned in the GeoTIFF driver docs though I don't believe it is specific about the interpretation of the resolution units. 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 Software Developer -- ---+-- 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 Software Developer ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_translate to TIFF -- color mismatch
Wow, that was fast. Thanks! I see that the fix is in 1.8.2... Do you know will this be a separate release or will the fixes be merged into 1.9.0? I ask because I see some recent discussions suggesting that 1.9.0 may be released by the end of the year. Thanks again for the quick response, Mladen Quoting Even Rouault even.roua...@mines-paris.org: Le jeudi 15 décembre 2011 18:08:26, mlade...@distributel.net a écrit : Hi all, I am trying to use gdal_translate to convert a .grd file (Northwood/VerticalMapper format) to TIFF. However, the colors in the translated file don't seem to be interpolated correctly based on the color scheme, for values near the end of the range. Here is a link to the original GRD, as well as the exported TIFFs (small download, only 200K): http://dl.dropbox.com/u/53500018/GdalQuestion.zip The GdalTranslateExport.tif file was produced with the command gdal_translate -b 1 -b 2 -b 3 C:\source\Ottawa_east_DEM30m_trim_copy.grd C:\temp\GdalTranslateExport.tif. Notice the top and bottom areas of the image: they are yellow. For comparison, take a look at the NativeExport.tif file in the same package. The top and bottom areas are orange, as they should be. The NativeExport.tif file was produced by Pitney Bowes MapInfo. Here is the color scheme for the GRD: http://www.smugmug.com/photos/i-bRZ879w/0/X3/i-bRZ879w-X3.jpg Notice that the thresholds are 70 for yellow, and 161 for orange. Well, the values at the top and bottom are in the high 130s/low 140s. So they should be interpolated close to orange, not yellow. Therefore the NativeExport.tif is correct, and the gdal_translate export is incorrect. I have also tried changing the thresholds, and as soon as I change the orange threshold to be below the maximum actual value in the GRD (let's say 144), the gdal_translate export shows the orange shades. Good analysis. There was indeed a bug in the way the end of the color palette is interpolated if the maximum z value is below a threshold of the color scheme. See http://trac.osgeo.org/gdal/ticket/4395 for the fix (a simple typo) So, it seems there might be a bug in gdal_translate where color values are not interpolated correctly for thresholds that fall outside of the range of actual values? Or is there a way to make gdal_translate behave the way that I want with some different command line options? Regards, Mladen Gavrilovic P.S. For convenience, here is a gallery with the TIFFs converted to JPEG (difference is still obvious), and the color scheme: http://smu.gs/vMEemF P.P.S Version of GDAL is 1.8.1 This message was sent using Distributel Webmail. Ce message a été envoyé à partir de la Messagerie Web Distributel. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev This message was sent using Distributel Webmail. Ce message a été envoyé à partir de la Messagerie Web Distributel. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] gdal_translate to TIFF -- color mismatch
Le jeudi 15 décembre 2011 23:32:14, mlade...@distributel.net a écrit : Wow, that was fast. Thanks! I see that the fix is in 1.8.2... Do you know will this be a separate release or will the fixes be merged into 1.9.0? I ask because I see some recent discussions suggesting that 1.9.0 may be released by the end of the year. The fix also went to trunk (1.9.0) of course (*). And 1.9.0 will likely be released before 1.8.2 indeed. (*) Trac has no way to specify multiple target milestones, so I tend to indicate the milestone corresponding to the stable branch if it has been applied there too. The general rule being that all fixes that go to the stable branch go to trunk too. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Slope in gdaldem
Here's what I've settled on. Let me know your thoughts. The algorithm now correctly calculates angle to lowest nearby point. If it's a local minima, it instead calculates angle to the nearest high point, which will be a negative number, so if the user wants to nullify these or use them in additional processing steps they can. How does adding code to gdal work? Best, Steve float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells const double radiansToDegrees = 180.0 / M_PI; GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData; double key; float pafLocalMin = afWin[0]; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMin) { pafLocalMin=afWin[k]; } if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } key = afWin[4] - pafLocalMin; if (key 0) { key = afWin[4] - pafLocalMax; } if (psData-slopeFormat == 1) return (float) (atan(sqrt(key) / (2*psData-scale)) * radiansToDegrees); else return (float) (100*(sqrt(key) / (2*psData-scale))); } ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Slope in gdaldem
float pafLocalMin = afWin[0]; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMin) { pafLocalMin=afWin[k]; } if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } key = afWin[4] - pafLocalMin; if (key 0) { key = afWin[4] - pafLocalMax; } Stephen 'key' can never be negative because pafLocalMin is computed from 'afWin' alone. The min val it can have is 0, when the minimum value is the central node itself. Joaquim ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Slope in gdaldem
Hmm. I don't follow? The revision would be this: float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells const double radiansToDegrees = 180.0 / M_PI; GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData; double key; float pafLocalMin = afWin[0]; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMin) { pafLocalMin=afWin[k]; } if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } key = afWin[4] - pafLocalMin; if (psData-slopeFormat == 1) return (float) (atan(sqrt(key) / (2*psData-scale)) * radiansToDegrees); else return (float) (100*(sqrt(key) / (2*psData-scale))); } But, clearly I'm not understanding something about the window. Steve On Thu, Dec 15, 2011 at 10:05 PM, Joaquim Luis jl...@ualg.pt wrote: float pafLocalMin = afWin[0]; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMin) { pafLocalMin=afWin[k]; } if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } key = afWin[4] - pafLocalMin; if (key 0) { key = afWin[4] - pafLocalMax; } Stephen 'key' can never be negative because pafLocalMin is computed from 'afWin' alone. The min val it can have is 0, when the minimum value is the central node itself. Joaquim ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Slope in gdaldem
On 16-12-2011 03:31, Stephen Mather wrote: Hmm. I don't follow? The revision would be this: float GDALSlopeHydroAlg (float* afWin, float fDstNoDataValue, void* pData) { // Hydrologic Slope is the max // local slope btw center cell and adjacent cells const double radiansToDegrees = 180.0 / M_PI; GDALSlopeAlgData* psData = (GDALSlopeAlgData*)pData; double key; float pafLocalMin = afWin[0]; float pafLocalMax = afWin[0]; for ( int k = 1; k 9; k++) { if (afWin[k] pafLocalMin) { pafLocalMin=afWin[k]; } if (afWin[k] pafLocalMax) { pafLocalMax=afWin[k]; } } key = afWin[4] - pafLocalMin; if (psData-slopeFormat == 1) return (float) (atan(sqrt(key) / (2*psData-scale)) * radiansToDegrees); else return (float) (100*(sqrt(key) / (2*psData-scale))); } But, clearly I'm not understanding something about the window. Steve, there are more things I do not understand either. What is 'psData-scale'? Is it the grid cell size? So, if 'pafLocalMin' is on one of the 4 corners the distance is not cell size (let's call it DX) but instead sqrt(DX*DX + DX*DX). And what if data is in geogs? Than DX is clearly different from DY and changes with latitude. You need to take that into account to compute the slope. Joaquim ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: Slope in gdaldem
Steve, there are more things I do not understand either. What is 'psData-scale'? Is it the grid cell size? So, if 'pafLocalMin' is on one of the 4 corners the distance is not cell size (let's call it DX) but instead sqrt(DX*DX + DX*DX). And what if data is in geogs? Than DX is clearly different from DY and changes with latitude. You need to take that into account to compute the slope. :) I don't plan to refactor this for non-planar coordinate systems. That is beyond my skills. AFAIK, none of the gdaldem code would be appropriate to run on lat lon-- which is probably something which would be good to address, since many global datasets default to that, but I can barely read c++, let along do that level of refactoring. Of course, if I'm wrong in how I'm reading the source for the existing gdaldem code, and it does deal with non-planar cases, I'll adapt that code for this use case. As to corners and sqrt(2), I'll take another look through tomorrow. In the likely case that you're right, I was sloppy in adapting the roughness index for this purpose and will have to take those distances into account. This will of course change the metric by which the lowest nearby point is determined, because of how distance will affect the factor of slope. And so, my 10 minute solution takes much longer in the end... . Best, Steve ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Runtime error on executing OGRRegisterAll() function.
Thanks. I've done things you have written, but it's still same. -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/Runtime-error-on-executing-OGRRegisterAll-function-tp4453903p7099963.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