Re: [gdal-dev] Support USRP

2013-11-18 Thread xavier lhomme
Hi

  What should be the return of gdalinfo after reading a TRANSH01.THF :
   I have added all IMG in subdatasets. But should I set a SpatialReference
and Extent for the main dataset (by merging all Extent in one ) ?

xav


2013/11/15 Even Rouault even.roua...@mines-paris.org

 Le jeudi 14 novembre 2013 10:35:35, xavier lhomme a écrit :
  Hi
 
I have a USRP product with several images organized with
 a THF file which describe all containing image and a directory for
 each
  image
  (with  .Gen, .SOU, .IMG and a .QAL ).
  The ASRP/USRP driver recognize .IMG file as explain in the web site.
 
  ASRP/USRP is related to ADRG (and source code is stored in the ADRG
  directory)
  the ADRG driver is able to recognize THF file
 
  Is it possible to improve ASRP/USRP driver to support THF files ?

 That is undoubtedly technically doable. The THF being a catalog, the best
 use
 of it I can see in an enhanced driver would be to report those images as
 SUBDATASETS of the THF dataset.

 
 
  Best regards
  xlhomme

 --
 Geospatial professional services
 http://even.rouault.free.fr/services.html

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Support USRP

2013-11-18 Thread Even Rouault
Selon xavier lhomme lhomme.xav...@gmail.com:

 Hi

   What should be the return of gdalinfo after reading a TRANSH01.THF :
I have added all IMG in subdatasets. But should I set a SpatialReference
 and Extent for the main dataset (by merging all Extent in one ) ?

Datasets that only report subdatasets generally don't define any spatial
reference or extent. Basically they keep with the defaults of the GDALDataset
constructor.


 xav


 2013/11/15 Even Rouault even.roua...@mines-paris.org

  Le jeudi 14 novembre 2013 10:35:35, xavier lhomme a écrit :
   Hi
  
 I have a USRP product with several images organized with
  a THF file which describe all containing image and a directory for
  each
   image
   (with  .Gen, .SOU, .IMG and a .QAL ).
   The ASRP/USRP driver recognize .IMG file as explain in the web site.
  
   ASRP/USRP is related to ADRG (and source code is stored in the ADRG
   directory)
   the ADRG driver is able to recognize THF file
  
   Is it possible to improve ASRP/USRP driver to support THF files ?
 
  That is undoubtedly technically doable. The THF being a catalog, the best
  use
  of it I can see in an enhanced driver would be to report those images as
  SUBDATASETS of the THF dataset.
 
  
  
   Best regards
   xlhomme
 
  --
  Geospatial professional services
  http://even.rouault.free.fr/services.html
 



___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] Support USRP

2013-11-18 Thread xavier lhomme
Hi
 I've submitted a ticket #5297 and attached a new version of the usrp/asrp
driver.

Best Regards


2013/11/18 Even Rouault even.roua...@mines-paris.org

 Selon xavier lhomme lhomme.xav...@gmail.com:

  Hi
 
What should be the return of gdalinfo after reading a TRANSH01.THF :
 I have added all IMG in subdatasets. But should I set a
 SpatialReference
  and Extent for the main dataset (by merging all Extent in one ) ?

 Datasets that only report subdatasets generally don't define any spatial
 reference or extent. Basically they keep with the defaults of the
 GDALDataset
 constructor.

 
  xav
 
 
  2013/11/15 Even Rouault even.roua...@mines-paris.org
 
   Le jeudi 14 novembre 2013 10:35:35, xavier lhomme a écrit :
Hi
   
  I have a USRP product with several images organized with
   a THF file which describe all containing image and a directory for
   each
image
(with  .Gen, .SOU, .IMG and a .QAL ).
The ASRP/USRP driver recognize .IMG file as explain in the web site.
   
ASRP/USRP is related to ADRG (and source code is stored in the ADRG
directory)
the ADRG driver is able to recognize THF file
   
Is it possible to improve ASRP/USRP driver to support THF files ?
  
   That is undoubtedly technically doable. The THF being a catalog, the
 best
   use
   of it I can see in an enhanced driver would be to report those images
 as
   SUBDATASETS of the THF dataset.
  
   
   
Best regards
xlhomme
  
   --
   Geospatial professional services
   http://even.rouault.free.fr/services.html
  
 



___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] [PATCH v2] Support Mercator_2SP in GeoTIFF

2013-11-18 Thread Antti Castrén
Hi Trent,

The chart  opens in ArcMap (ArcGIS 10.0) well, and it is in right
location (Seattle).

Relevant properties of the file as seen by ArcGIS:
Cell Size: 240.003787, 240.003787
Extent
Top:4143530.20898
Left: -9278434.53415
Right:  -9165872.75807
Bottom:  3984167.69444
Spatial Reference: Global Mercator
Linear Unit: Meter
false_easting: 0
false_northing: 0
central_meridian: 0
standard_parallel_1:  47.667
Datum: D_North_American_1983


I opened the file in FME Data Inspector (FME Desktop 2013 SP1) also.
It did not read the latitude of true scale from the
ProjStdParallel1GeoKey as you can see from the following parameters.
Therefore it places the chart to NE Georgia/NW South Carolina as if
the latitude of true scale was the Equator.

FME Parameters
Coordinate System Parameters
CS_NAME: _FME_0
DESC_NM: Global Mercator
DT_NAME: NAD83
PROJ: MRCAT
UNIT: METER

Datum Parameters
DESC_NM: NAD 1983, Alaska, Canada,Continental US,Mexico,Central America
ELLIPSOID: GRS1980
SOURCE: US Defense Mapping Agency, TR-8350.2-B, August 1993
USE: NAD83

Ellipsoid Parameters
DESC_NM: Geodetic Reference System of 1980
E_RAD: 6378137
P_RAD: 6356752.3141403478
SOURCE: Stem, L.E., Jan 1989, State Plane Coordinate System of 1983

OGC WKT Description
PROJCS[Global Mercator,
GEOGCS[NAD83,
DATUM[North_American_Datum_1983,
SPHEROID[Geodetic Reference System of
1980,6378137,298.2572221008916,
AUTHORITY[EPSG,7019]],
AUTHORITY[EPSG,6269]],
PRIMEM[Greenwich,0],
UNIT[degree,0.0174532925199433],
AUTHORITY[EPSG,4269]],
PROJECTION[Mercator_1SP],
PARAMETER[central_meridian,0],
PARAMETER[scale_factor,1],
PARAMETER[false_easting,0],
PARAMETER[false_northing,0],
UNIT[METER,1]]

Esri WKT Description
PROJCS[Global_Mercator,GEOGCS[GCS_North_American_1983,DATUM[D_North_American_1983,SPHEROID[Geodetic_Reference_System_of_1980,6378137,298.2572221008916]],PRIMEM[Greenwich,0],UNIT[Degree,0.017453292519943295]],PROJECTION[Mercator],PARAMETER[central_meridian,0],PARAMETER[false_easting,0],PARAMETER[false_northing,0],UNIT[Meter,1],PARAMETER[standard_parallel_1,0.0]]
/FME Parameters

Intergraph's Geomedia Professional recognises the file as GeoTIFF, but
the chart appears to be in the same false place as with FME.

If you want, I could update the file's georeferencing by FME and
Geomedia so that location is correct. Then I could send listgeo
information or even put the actual files available if desired.

The goal is, of course, to get various software vendors to use the
GeoTIFF-information the same way. Hopefully the common way would also
be logically correct, but I guess that's too much to ask.

Best regards,

  Antti


Here is the listgeo info from the original file for those who are interested:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
  ModelTiepointTag (2,3):
 000
 -9278434.53  4143530.21   0
  ModelPixelScaleTag (1,3):
 240.003787   240.003787   0
  End_Of_Tags.
   Keyed_Information:
  GTModelTypeGeoKey (Short,1): ModelTypeProjected
  GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
  GTCitationGeoKey (Ascii,16): Global Mercator
  GeographicTypeGeoKey (Short,1): GCS_NAD83
  GeogCitationGeoKey (Ascii,6): NAD83
  GeogAngularUnitsGeoKey (Short,1): Angular_Degree
  GeogSemiMajorAxisGeoKey (Double,1): 6378137
  GeogInvFlatteningGeoKey (Double,1): 298.257222
  Unknown-2062 (Double,3): 000
  ProjectedCSTypeGeoKey (Short,1): User-Defined
  ProjectionGeoKey (Short,1): User-Defined
  ProjCoordTransGeoKey (Short,1): CT_Mercator
  ProjLinearUnitsGeoKey (Short,1): Linear_Meter
  ProjStdParallel1GeoKey (Double,1): 47.667
  ProjNatOriginLongGeoKey (Double,1): 0
  ProjNatOriginLatGeoKey (Double,1): 0
  ProjFalseEastingGeoKey (Double,1): 0
  ProjFalseNorthingGeoKey (Double,1): 0
  ProjScaleAtNatOriginGeoKey (Double,1): 1
  End_Of_Keys.
   End_Of_Geotiff.
Projection Method: CT_Mercator
   ProjNatOriginLatGeoKey: 0.00 (  0d 0' 0.00N)
   ProjNatOriginLongGeoKey: 0.00 (  0d 0' 0.00E)
   ProjScaleAtNatOriginGeoKey: 1.00
   ProjFalseEastingGeoKey: 0.00 m
   ProjFalseNorthingGeoKey: 0.00 m
GCS: 4269/NAD83
Datum: 6269/North American Datum 1983
Ellipsoid: 7019/GRS 1980 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.00/  0d 0' 0.00E)
Projection Linear Units: 9001/metre (1.00m)
Corner Coordinates:
Upper Left(-9278434.534, 4143530.209)
Lower Left(-9278434.534, 3984167.694)
Upper Right   (-9165872.758, 4143530.209)
Lower Right   (-9165872.758, 3984167.694)
Center(-9222153.646, 4063848.952)


2013/11/14 Trent Piepho tpie...@gmail.com:
 It looks like arcgis does not support Mercator_1SP in geotiff files.
 It just ignores the scale factor and uses 1.0.  Which means in GDAL
 and libgeotiff as they are 

[gdal-dev] Patch: Improvements to the rasdaman driver

2013-11-18 Thread Fabian Schindler
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Devs,

I implemented a couple of improvements of the GDAL rasdaman driver.
Ticket with a detailed description of the changes and patch are
available here: http://trac.osgeo.org/gdal/ticket/5298


Of course, I appreciate any feedback and I hope the patch will be applied!

Cheers,
Fabian
-BEGIN PGP SIGNATURE-
Version: GnuPG v1.4.14 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/

iQEcBAEBAgAGBQJSijTpAAoJEJoaok7HgajKsXsH/0bjVpeaH1CSanzEvYQrAqpt
fGbv5mfNPoJ/CIh+eiA12UWVqOfkEUhsLsOHqKdPyWxTFUHt2++U4mBKyk0LbtsX
ihJZFoZ0VGnrwqvoHgbha44f7QIX5iDTsuVuMIN6Q+8cNOX1vckyGMDMzSmdSGDR
2mmm80nAtOxdVF5OwMaOXBCF2DraVMPjk0xd9sntSGJ1EgQgfxUP2p5w7ItTP+UG
pAw0xK+D+0kpoj4vkyq2qWWTkC56O1GLUoZP8s1WkdhxMWRCtkrKKevDpdStg2NU
kPtLu0/fc2IqU9d/lL4PJhlS7kMjJ5wz4XWuJ0whwlAx27MNWH6L5P1bdaSCyWs=
=MOCq
-END PGP SIGNATURE-
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev


[gdal-dev] ogr2og / forcing a +-180 system

2013-11-18 Thread Robb K. Wright


I need to use ogr2ogr (or any other command line) to convert a shapefile 
with longitudes spanning from -220 to -60, forcing it into a +-180 
scheme, so the -220 coordinate would come out as +140 and the -60 stays 
as -60.  I don't need to worry about the polys that actually cross the 
180 line.


Any thoughts?  I was thinking that running it through EPSG:4326 would 
reshape it (as it does in Esri-land), but no luck.


Robb



___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] ogr2og / forcing a +-180 system

2013-11-18 Thread Frank Warmerdam
Robb,

Have you tried the -wrapdateline switch?  I do not believe geometries that
cross the dateline will be handled ideally (they aren't split as one might
hope), but if you don't have them then the rest of the geometries should be
wrapped.  I see there is even now a -datelineoffset switch if someone
needed 0 to 360.

Best regards,
Frank


On Mon, Nov 18, 2013 at 8:36 AM, Robb K. Wright robbkwri...@gmail.comwrote:


 I need to use ogr2ogr (or any other command line) to convert a shapefile
 with longitudes spanning from -220 to -60, forcing it into a +-180 scheme,
 so the -220 coordinate would come out as +140 and the -60 stays as -60.  I
 don't need to worry about the polys that actually cross the 180 line.

 Any thoughts?  I was thinking that running it through EPSG:4326 would
 reshape it (as it does in Esri-land), but no luck.

 Robb



 ___
 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

[gdal-dev] Generated KML incorrectly formated and image reflected in X axis

2013-11-18 Thread JSz
Good evening,

I have generated an KMZ file using :
gdal_translate.exe -of KMLSUPEROVERLAY EM_RESULTS.tiff EM_RESULTS.kmz -co
FORMAT=JPEG

Opening the KMZ on its own with GoogleEarth doesnt display anything, if I
extract the tiles directly out the KML and extract to Desktop and drill down
I can open the KML / images.  However the image when overlayed in GE seems
to be reflected in the X axis.
 
http://osgeo-org.1560.x6.nabble.com/file/n5089715/Capture.png 


gdal_grid -a invdist:smoothing=0:radius1=0.0004:radius2=0.0004 -of GTiff
-zField hcp_cond -ot Float64 -l data EMSCANNING.vrt EM_RESULTS.tiff


OGRVRTDataSource
OGRVRTLayer name=data
SrcDataSourcedata.csv/SrcDataSource
GeometryTypewkbPoint/GeometryType
LayerSRSWGS84/LayerSRS
GeometryField separator=  encoding=PointFromColumns x=lon y=lat
z=hcp_cond/
/OGRVRTLayer
/OGRVRTDataSource




Any pointers or suggestions on resolving both the projection and KMZ
structure?




--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Generated-KML-incorrectly-formated-and-image-reflected-in-X-axis-tp5089715.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] ogr2og / forcing a +-180 system

2013-11-18 Thread Robb K. Wright


I haven't been able to get either -wrapdateline or -datelineoffset to 
alter the coords.  As far as I can figure out, the -wrapdateline option 
only operates on features that actually intersect the 180 line, so using 
it doesn't alter my coords since mine are only on either side.  I'm just 
trying to shift any polys -180 to be have long +360.  Here's a sample 
shape if anybody wants to take a swing:  
http://branewright.org/test/test_file.zip


Robb

On 11/18/13 12:50 PM, Frank Warmerdam wrote:

Robb,

Have you tried the -wrapdateline switch?  I do not believe geometries 
that cross the dateline will be handled ideally (they aren't split as 
one might hope), but if you don't have them then the rest of the 
geometries should be wrapped.  I see there is even now a 
-datelineoffset switch if someone needed 0 to 360.


Best regards,
Frank


On Mon, Nov 18, 2013 at 8:36 AM, Robb K. Wright robbkwri...@gmail.com 
mailto:robbkwri...@gmail.com wrote:



I need to use ogr2ogr (or any other command line) to convert a
shapefile with longitudes spanning from -220 to -60, forcing it
into a +-180 scheme, so the -220 coordinate would come out as +140
and the -60 stays as -60.  I don't need to worry about the polys
that actually cross the 180 line.

Any thoughts?  I was thinking that running it through EPSG:4326
would reshape it (as it does in Esri-land), but no luck.

Robb



___
gdal-dev mailing list
gdal-dev@lists.osgeo.org mailto: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 mailto:warmer...@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam 
http://pobox.com/%7Ewarmerdam

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] Generated KML incorrectly formated and image reflected in X axis

2013-11-18 Thread Even Rouault
Le lundi 18 novembre 2013 19:14:15, JSz a écrit :
 Good evening,
 
 I have generated an KMZ file using :
 gdal_translate.exe -of KMLSUPEROVERLAY EM_RESULTS.tiff EM_RESULTS.kmz -co
 FORMAT=JPEG
 
 Opening the KMZ on its own with GoogleEarth doesnt display anything, if I
 extract the tiles directly out the KML and extract to Desktop and drill
 down I can open the KML / images.  However the image when overlayed in GE
 seems to be reflected in the X axis.
 
 http://osgeo-org.1560.x6.nabble.com/file/n5089715/Capture.png
 
 
 gdal_grid -a invdist:smoothing=0:radius1=0.0004:radius2=0.0004 -of GTiff
 -zField hcp_cond -ot Float64 -l data EMSCANNING.vrt EM_RESULTS.tiff

gdal_grid generates by default images with bottom-left corner as origin, 
whereas the usual convention is top-left. So this is perhaps your issue. You 
could try to do : gdalwarp EM_RESULTS.tiff EM_RESULTS_warped.tiff and 
gdal_translate this intermediary file to KMLSUPEROVERLAY.

 
 
 OGRVRTDataSource
 OGRVRTLayer name=data
 SrcDataSourcedata.csv/SrcDataSource
 GeometryTypewkbPoint/GeometryType
 LayerSRSWGS84/LayerSRS
 GeometryField separator=  encoding=PointFromColumns x=lon y=lat
 z=hcp_cond/
 /OGRVRTLayer
 /OGRVRTDataSource
 
 
 
 
 Any pointers or suggestions on resolving both the projection and KMZ
 structure?
 
 
 
 
 --
 View this message in context:
 http://osgeo-org.1560.x6.nabble.com/Generated-KML-incorrectly-formated-and
 -image-reflected-in-X-axis-tp5089715.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

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] building gdal on windows -- how to hide internal libtiff symbols

2013-11-18 Thread kyle.sletmoe
Even,

Oh, OK I was not aware of this. Thanks for the help.

Regards,
Kyle



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/building-gdal-on-windows-how-to-hide-internal-libtiff-symbols-tp5089510p5089742.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] Generated KML incorrectly formated and image reflected in X axis

2013-11-18 Thread JSz
Perfect, that solved the projection issue.

Have you any ideas why the KMZ file isn't able to open in google earth but
if I remove the individual parts ( 0.kml and 0.png ) and open all is ok?

Alternativly is there a way to just generate the KML file directly rather
than with :

gdal_translate.exe -of KMLSUPEROVERLAY MY_TIFF_File.tiff Output.kmz -co
FORMAT=PNG

J



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/Generated-KML-incorrectly-formated-and-image-reflected-in-X-axis-tp5089715p5089752.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] ogr2og / forcing a +-180 system

2013-11-18 Thread Even Rouault
Le lundi 18 novembre 2013 20:29:07, Robb K. Wright a écrit :
 I haven't been able to get either -wrapdateline or -datelineoffset to
 alter the coords.  As far as I can figure out, the -wrapdateline option
 only operates on features that actually intersect the 180 line, so using
 it doesn't alter my coords since mine are only on either side.  I'm just
 trying to shift any polys -180 to be have long +360.  Here's a sample
 shape if anybody wants to take a swing:
 http://branewright.org/test/test_file.zip

Robb,

I can think of 2 ways :

- doing an intermediate reprojection to EPSG:3857, and then reprojecting the 
result to EPSG:4326.

ogr2ogr tmp.shp test5.shp -t_srs EPSG:3857
ogr2ogr test5_shifted.shp tmp.shp -t_srs EPSG:4326

- with GDAL 1.10 and Spatialite 4.0, using the spatialite ST_Shift_Longitude() 
function :

ogr2ogr test5_shifted.shp test5.shp -dialect sqlite \
   -sql select ST_Shift_Longitude(geometry) as GEOMETRY, * from test5 

Even

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] set authenticate set delivery off

2013-11-18 Thread Kyle Sletmoe
set authenticate set delivery off

-- 
*Information contained herein is subject to the Code of Federal Regulations 
Chapter 22 International Traffic in Arms Regulations. This data may not be 
resold, diverted, transferred, transshipped, made available to a foreign 
national within the United States, or otherwise disposed of in any other 
country outside of its intended destination, either in original form or 
after being incorporated through an intermediate process into other data 
without the prior written approval of the US Department of State.  **Penalties 
for violation include bans on defense and military work, fines and 
imprisonment.*
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] How to implement tile read / write to gdal supported format file?

2013-11-18 Thread Luke
You could do something like the following:
(code modified from the gdal_calculations library -
http://code.google.com/p/gdal-calculations)

from osgeo import gdal

class Block(object):

def __init__(self, dataset, band, xoff, yoff, xsize, ysize):
self.xoff = xoff
self.yoff = yoff
self.xsize = xsize
self.ysize = ysize
if not band: self.data = dataset.ReadAsArray(xoff, yoff, xsize,
ysize)
else: self.data = dataset.GetRasterBand(band).ReadAsArray(xoff,
yoff, xsize, ysize)

def read_blocks_as_array(dataset, band=0, nblocks=1):
'''Read GDAL Datasets or individual Bands block by block'''

ncols=dataset.RasterXSize()
nrows=dataset.RasterYSize()
xblock,yblock=dataset.GetRasterBand(1).GetBlockSize()

if xblock==ncols:
yblock*=nblocks
else:
xblock*=nblocks

for yoff in xrange(0, nrows, yblock):

if yoff + yblock  nrows:
ysize = yblock
else:
ysize  = nrows - yoff

for xoff in xrange(0, ncols, xblock):
if xoff + xblock  ncols:
xsize  = xblock
else:
xsize = ncols - xoff

yield Block(dataset, band, xoff, yoff, xsize, ysize)


ds1=gdal.Open(some_raster)
ds2=gdal.Open(another_raster)
out_driver=gdal.GetDriverByName('GTIFF')
out_ds=out_driver.Create (out_raster,ds1.RasterXSize(),ds1.RasterYSize(),
   1,ds1.GetRasterBand(1).DataType,
   ['BIGTIFF=IF_SAFER'])

out_ds.SetGeoTransform(ds1.GetGeoTransform())
out_ds.SetProjection(ds1.GetProjection())
out_band=out_ds.GetRasterBand(1)

for block in read_blocks_as_array(ds1,band=1,nblocks=10):
block2=Block(ds2, 1, block.xoff, block.yoff, block.xsize, block.ysize)
ndarray1=block.data
ndarray2=block2.data
result = ndarray1+ndarray2
out_band.WriteArray(result, block.xoff, block.yoff)



--
View this message in context: 
http://osgeo-org.1560.x6.nabble.com/gdal-dev-How-to-implement-tile-read-write-to-gdal-supported-format-file-tp5089581p5089803.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