[gdal-dev] Building GDAL 1.7.3 Python bindings
Greetings, I am looking for some suggestions on how to build the python bindings for the latest gdal sources (1.7.3) for Windows. Environment: Windows 7, Visual Studio 2010, Python 2.5 Success: building gdal itself Failure: python module + bindings I keep hitting dead ends getting Python set up. I am trying to include details here because I searched and saw many "which version are you using" messages... sorry for the verbosity. I am not encountering ERROR messages, just trying to complete the build process and get up to date python bindings. I set up nmake.opt including PYDIR and PY_INST_DIR. I can't see anywhere those are USED but I followed docs and set them appropriately. (C:\Python25) And indeed, nothing pythonic appears to get built or installed. In desperation I tried just copying the python stuff into site-packages but of course that does not work, since there is no _gdal.dll built and thus Python cannot find it. I searched for hints on how to get there from here, but I only found a few other people with similar problems. I was able to build gdal with Visual Studio's "Build" command. I see references to building from the command line, but I don't have a clue how to do that in Windows. Something called nmake is referenced. I have not marched down that road yet. I got source code from http://trac.osgeo.org/gdal/wiki/DownloadSource I looked at I looked at http://pypi.python.org/pypi/GDAL/ which has lovely and usable Linux instructions but I am currently trapped in a Windows world at my otherwise wonderful day job. The text is duplicated in the source code under swig/python/Readme.txt The link on that page suggesting that I download binaries jumps to the old 1.5 version, but I can find this folder http://download.osgeo.org/gdal/win32/1.7 which seems to have only a couple plugins. As a Linux person it seems very odd to drill down 3 levels into source code only to find something like "Windows: go get the binaries somewhere else". Where do Windows binaries come from if they don't get built ever?? So mysterious. I suspect I need to do something with swig, which I downloaded and installed, but the "Windows" notes basically say "you should never have to do this". I don't know what I should never have to do, so that is as far as I have gotten. I can't use FWTools as they are built against an old version of Python (2.3) and I am trying to work in an ESRI world which means I am using 2.5 and 2.6 and further I really would like to be able to build support for file geodatabases to make this gdal python thing really usable in my work. I am sorry this email is long, but I have had plenty of time to work on building GDAL today while waiting for fragile and slow ESRI tools to run, really motivating me to try to use GDAL + Python. (Did the geoprocess crash or is it just sloo... one never knows with ESRI Model Builder. OH! The process just FINISHED. Unreal. I had time to download and build gdal itself a couple times while waiting. Now I can go back to real work.) Cheers -- Brian Wilson Corvallis Oregon ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] .asc files during open
Best way is probably to add a driver option (if it doesn't already have one for that). Ray On 12/3/2010 3:01 PM, Todd Smith wrote: When opening .asc files, they are scanned for any values that may be floating point so the type of data (int or float) can be set properly. For the case I'm dealing with (2800 files @ 50mb each) it means that each file takes about 1.25 seconds to open (because the data is all integer data). It means it takes nearly an hour for all the files to be opened. Short of changing the code to not scan for a '.' and just assume that every file is a floating point file, can anyone think of a good way to speed this up? Thanks, Todd ___ 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] .asc files during open
When opening .asc files, they are scanned for any values that may be floating point so the type of data (int or float) can be set properly. For the case I'm dealing with (2800 files @ 50mb each) it means that each file takes about 1.25 seconds to open (because the data is all integer data). It means it takes nearly an hour for all the files to be opened. Short of changing the code to not scan for a '.' and just assume that every file is a floating point file, can anyone think of a good way to speed this up? Thanks, Todd ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] several No data value
Ludovic Mercier wrote: Hi, i work with a type of data wich can have : no data value and other information like "low and hight instrument saturation". For many algorithm this value (instrument saturation ) is like "no data value". But in gdal doc only one data value is managed. Have you a solution for this case ? A practical exemple is : data is REAL typically [-2E-5, 2E-5] no data value is -1004 instrument saturation is -1001 thanks by advance for yours suggestions. Ludovic, GDAL only supports single nodata values as part of the data model. You can add arbitrary metadata describing other nodata values, but generally GDAL itself and other client software won't realize the significance of these other values. sorry! 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
Re: [gdal-dev] GDAL/OGR Python Question
Brian Walawender wrote: Hello, I am looking for some assistance in optimizing a section of code for faster performance. Here is my problem, I have GeoTiff that contains 1 km x 1 km population data over an area roughly the size of the continental US (7020 x 3000). I am trying to calculate the population within a polygon. I am doing this by finding the extent of the polygon and reading in that section of the GeoTiff using the ReadAsArray function. At this point I can quickly calculate the population within the extent by using the numpy sum. However, the only way I can figure to sum the points within the polygon is to iterate over the array checking each point using ogr. If the polygon is very large, this can take an extended period of time. Is there a faster way to calculate the population within a large polygon? My code and some sample output is below. Brian, I think you would be better off rasterizing your mask polygon and operating on the masked raster. It might be sufficient to set all cells outside the region to 0 in your target raster if you are just summing. Or you might generate the mask separately and just check it instead of doing an expensive point in polygon test. The test script for the rasterize function is available at the following url and might provide a hint of how to use it. http://svn.osgeo.org/gdal/trunk/autotest/alg/rasterize.py You can also find more about the rasterize api from it's C++ docs at: http://www.gdal.org/gdal__alg_8h.html#50caf4bc34703f0bcf515ecbe5061a0a Good luck, -- ---+-- 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
Re: [gdal-dev] GDAL contouring problem
FYI, I've commited a fix in GDAL trunk so that 0 as a floating point value can
be used as a valid value for -fl
Le mercredi 01 décembre 2010 05:44:09, Ole Nielsen a écrit :
> Thanks for explaining the reason why the numeric value of 0.00 doesn't work
> as a fixed level in gdal_contour while the value of 0 does. I understand
> now the reason for this.
> However I would have thought one could replace the test for 'zeroness' in
> (atof(argv[i+1]) != 0 || EQUAL(argv[i+1], "0"))
> with a test for whether argv]i+1] is not empty and can be converted to a
> numeric value. This would be a duck-typing approach and more robust in my
> humble opinion.
>
> Best regards
> Ole Nielsen
>
>
> -Original Message-
> From: gdal-dev-boun...@lists.osgeo.org
> [mailto:gdal-dev-boun...@lists.osgeo.org] On Behalf Of
> gdal-dev-requ...@lists.osgeo.org Sent: Monday, 29 November 2010 8:41 PM
> To: gdal-dev@lists.osgeo.org
> Subject: gdal-dev Digest, Vol 78, Issue 64
>
> Send gdal-dev mailing list submissions to
> gdal-dev@lists.osgeo.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
> or, via email, send a message with subject or body 'help' to
> gdal-dev-requ...@lists.osgeo.org
>
> You can reach the person managing the list at
> gdal-dev-ow...@lists.osgeo.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of gdal-dev digest..."
>
>
> Today's Topics:
>
>1. Re: GDAL contouring problem (Frank Warmerdam)
>2. Re: GDAL contouring problem (Francis Markham)
>3. Re: Delete a sqlite database (Ludovic Granjon)
>4. Re: Error building GDAL with GEOS on Solaris (Namrata)
>5. Making non-rectangular image edges transparent
> (Just van den Broecke)
>6. geotiff conflict (Julien Malik)
>7. Re: RFC 31 - OGR 64bit Support (David Burken)
>
>
> --
>
> Message: 1
> Date: Sun, 28 Nov 2010 23:24:56 -0500
> From: Frank Warmerdam
> Subject: Re: [gdal-dev] GDAL contouring problem
> To: Ole Nielsen
> Cc: "gdal-dev@lists.osgeo.org" ,
> "adeleb...@me.com"
> Message-ID: <4cf32b18.2030...@pobox.com>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> Ole Nielsen wrote:
> > I am posting this again without test files attached due the 50K size
> > limit. Anyone who wants to see them, please contact me directly.
> >
> >
> >
> > --
> >
> > We have observed interesting anomaly with gdal_contour when one of the
> > fixed levels is zero.
> >
> > If the zero contour is expressed as 0.0 (or indeed 0.) gdal_contour
> > replies with the standard Usage message (see below). If the zero contour
> > is expressed as the integer 0 (or 0.0001) it works and produces the
> > expected contours.
>
> Ole,
>
> The gdal_contour program has some rather hokey heuristics to try and
> recognise the end of the list of contour levels. The code looks like:
>
> else if( EQUAL(argv[i],"-fl") && i < argc-1 )
> {
> while( i < argc-1
> && nFixedLevelCount
> <
> (int)(sizeof(adfFixedLevels)/sizeof(double)) && (atof(argv[i+1]) != 0 ||
> EQUAL(argv[i+1],"0")) && !EQUAL(argv[i+1], "-3d"))
> adfFixedLevels[nFixedLevelCount++] = atof(argv[++i]);
> }
>
> So basically, it assumes everything is a "level" until something that
> has a numeric value of 0 is encounter that isn't the specific string
> "0". So your analysis of the problem is essentially correct, but it is
> more or less intentional as we try to support a list of levels with no
> explicit "end of list" marker in the arguments to the command.
>
> 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
>
>
>
> --
>
> Message: 2
> Date: Mon, 29 Nov 2010 16:48:56 +1100
> From: Francis Markham
> Subject: Re: [gdal-dev] GDAL contouring problem
> To: Frank Warmerdam
> Cc: "gdal-dev@lists.osgeo.org" , Ole Nielsen
> , "adeleb...@me.com"
> Message-ID:
>
> Content-Type: text/plain; charset=ISO-8859-1
>
> Is there any reason not to include a heuristic that uses strtod() or
> similar to test for the zeroness of the last argument?
>
> -Francis Markham
>
> On 29 November 2010 15:24, Frank Warmerdam wrote:
> > Ole Nielsen wrote:
> >> I am posting this again without test files attached due the 50K size
> >> limit. Anyone who wants to see them, please contact me directly.
> >>
> >>
> >> --
> >>
> >> We have observed interesting anomaly with gdal_contour when one of the
> >> fixed levels is zero.
> >>
> >> If the zero contour
[gdal-dev] GDAL/OGR Python Question
Hello,
I am looking for some assistance in optimizing a section of code for faster
performance. Here is my problem, I have GeoTiff that contains 1 km x 1 km
population data over an area roughly the size of the continental US (7020 x
3000). I am trying to calculate the population within a polygon. I am
doing this by finding the extent of the polygon and reading in that section
of the GeoTiff using the ReadAsArray function. At this point I can quickly
calculate the population within the extent by using the numpy sum. However,
the only way I can figure to sum the points within the polygon is to iterate
over the array checking each point using ogr. If the polygon is very large,
this can take an extended period of time. Is there a faster way to
calculate the population within a large polygon? My code and some sample
output is below.
Thanks.
bw
#findPopulationOGR.py - read a polygon from the database
# and calculate the population within it
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo.gdalconst import *
from osgeo.gdal_array import *
import sys, time
import numpy as np
populationFile = "pop4269.tif"
band = 1
# Zone based
#segment_id = 147231
# Polygon based
segment_id = 147810
sql = "SELECT AsBinary(ST_MemUnion(COALESCE (ppoly.geom, pl.geom))) AS
wkb_geometry "
sql = sql + "FROM iris.product_segment ps "
sql = sql + "JOIN iris.product_segment_to_product_location pspl ON ps.id =
pspl.product_segment_id "
sql = sql + "JOIN iris.product_location pl ON pspl.product_location_id =
pl.id "
sql = sql + "JOIN iris.product_location_type plt ON
pl.product_location_type_id = plt.id "
sql = sql + "LEFT JOIN iris.product_polygon ppoly ON ps.id =
ppoly.product_segment_id "
sql = sql + " WHERE ps.id=" + `segment_id`
conn = ogr.Open("PG: host='localhost' dbname=mydb' ")
layer = conn.ExecuteSQL(sql)
feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
#create memory layer for polygon
srWGS84 = osr.SpatialReference()
srWGS84.ImportFromProj4('+proj=longlat +datum=WGS84')
polyds = ogr.GetDriverByName('Memory').CreateDataSource('')
polylyr = polyds.CreateLayer('poly',geom_type=ogr.wkbPolygon, srs=srWGS84)
feat = ogr.Feature(polylyr.GetLayerDefn())
feat.SetGeometryDirectly(geom)
polylyr.CreateFeature(feat)
#check polygon layer extent
extent = polylyr.GetExtent()
print 'Polygon Bounding Box'
print 'UL: ', extent[0], extent[3]
print 'LR: ', extent[1], extent[2]
XUL = extent[0]
YUL = extent[3]
XLR = extent[1]
YLR = extent[2]
gdal.AllRegister()
ds = gdal.Open( populationFile, GA_ReadOnly )
if ds is None:
print('Cannot open %s' % populationFile)
sys.exit(1)
cols = ds.RasterXSize
rows = ds.RasterYSize
proj = ds.GetProjection()
geomatrix=ds.GetGeoTransform()
(success, inv_geometrix) = gdal.InvGeoTransform(geomatrix)
xUL = int(inv_geometrix[0] + inv_geometrix[1] * XUL + inv_geometrix[2] *
YUL)
yUL = int(inv_geometrix[3] + inv_geometrix[4] * XUL + inv_geometrix[5] *
YUL)
xLR = int(inv_geometrix[0] + inv_geometrix[1] * XLR + inv_geometrix[2] *
YLR)
yLR = int(inv_geometrix[3] + inv_geometrix[4] * XLR + inv_geometrix[5] *
YLR)
xdif = xLR - xUL
ydif = yUL - yLR
#print cols, rows, proj, xUL, yUL, xLR, yLR, xdif, ydif
data = ds.GetRasterBand(band).ReadAsArray(xUL, yLR, xdif, ydif)
extentPop = np.sum(data)
print "Print total population in extent = %#10d" % (extentPop )
startTime = time.time()
for i in range(data.shape[0]):
for j in range(data.shape[1]):
X = i + xUL
Y = j + yLR
lon = geomatrix[0] + X*geomatrix[1] + Y*geomatrix[2]
lat = geomatrix[3] + X*geomatrix[4] + Y*geomatrix[5]
qpt = ogr.CreateGeometryFromWkt("POINT(%s %s)" % (lon, lat))
if geom.Contains(qpt) == False:
data[i,j] = 0
polyPop = np.sum(data)
endTime = time.time()
print "Print total population in polygon = %#10d" % (polyPop )
print "The calc took " + str(endTime-startTime) + " seconds"
ds = None
conn = None
polylyr = None
$ python findPopulationOGR.py
Polygon Bounding Box
UL: -108.891024 49.22
LR: -104.041596 46.540403
Print total population in extent = 101602
Print total population in polygon = 21052
The calc took 604.560879946 seconds
--
Brian Walawender
Technique Development Meteorologist
Scientific Services Division - Central Region Headquarters
816-268-3114 - Office
816-805-6497 - Cell
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdal_translate and modis ocean HDFs
Hi everyone, I am trying to import some L3 MODIS satellite images into GRASS 6.4 by following the instructions from: http://grass.osgeo.org/wiki/MODIS#Import I've tried method 1 but I got some error messages. gdal_translate -a_srs "+init=epsg:4326" -a_nodata 65535 -a_ullr -180 90 180 -90 -co "COMPRESS=PACKBITS" HDF4_SDS:UNKNOWN:"A20100252010032.L3m_8D_CHL_chlor_a_9km":0 prueba.tif Input file size is 4320, 2160 0ERROR 1: SDreaddata() failed for block. ERROR 1: IReadBlock failed at X offset 0, Y offset 0 ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0 The error is the same with different files. My system is Ubuntu 9.10 64bits. gdal_translate --version GDAL 1.7.2, released 2010/04/23 What can I do? Thanks, Gustavo Martínez ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
