[gdal-dev] Question about compiling C++ gdal utility

[Roger André]

I'm trying to compile a gdal utility written in C++ (which is something
completely new to me), and am encountering an error.

color-relief.cpp:30:23: gdal_priv.h: No such file or directory
color-relief.cpp:31:23: stringtok.h: No such file or directory

The gdal_priv.h file is located in /export/include/gdal_priv.h, but I'm not
sure where the stringtok.h file is.

The command I gave to compile was this, sudo g++ color-relief.cpp -lgdal -o
color-relief.

Can one of you C++ geniuses give me a couple pointers about how to do this
the right way?

Thanks,

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

Re: [gdal-dev] HowTo open NetCDF SUBDATASET and write to new file

[Roger André]

Hi Steve,

Thanks very much for the code snippet.  I figured out the following via
experimentation.

- Open specific layer by using:
layer = gdal.Open('netCDF:annual_mean.nc:temp')

- Then get at the layer data data using the following:
netcdf_band = layer.GetRasterBand(1)

- Read in the data one row at a time:
for iY in range(y_size):
  netcdf_data = netcdf_band.ReadAsArray(0, iY, x_size, 1)

It works, but is a complete hack in the sense that I just plugged in stuff
from a bunch of different scripts and banged away at it until it gave me
what I needed.  Glad to see that it is similar to your code.  I'm curious
about the .reshape() method that you pass to ReadAsArray above.  What
does that do?

Thanks again for the help.

Roger
--

On Tue, Sep 16, 2008 at 12:05 AM, Steve Gaffigan [EMAIL PROTECTED]wrote:

 Hello.  Here's an example of reading an entire array for a single
 subdataset.

 [EMAIL PROTECTED] AEFF]$ gdalinfo latest.nc
 Driver: netCDF/Network Common Data Format
 ...
  SUBDATASET_25_NAME=NETCDF:latest.nc:TMP_GPML
  SUBDATASET_25_DESC=[49x11x78x141]
 air_temperature_at_constant_altitude_above_m
 ean_sea_level (8-bit integer)
 ...

 [EMAIL PROTECTED] AEFF]$ python
  from osgeo import gdal
  ds = gdal.Open('NETCDF:latest.nc:TMP_GPML')
  data = ds.ReadAsArray().reshape(49,11,78,141)

 Using a more generic call of the form below you can get at particular
 hyperslabs:

 
 data=ds.GetRasterBand(i).ReadAsArray(xoff=x0,yoff=y0,win_xsize=nx,win_ysize=ny)

 where, for the above example, i=[1,...,539].

 Steve


  Hi List,
 
  I'm would like to use the Python bindings to open a specific SUBDATASET
 in
  a
  NetCDF file and read the contents.  I read in
  http://www.gdal.org/gdal_datamodel.html, under the SUBDATASETS domain
  section, that I should be able to pass the _NAME parameter into
  GDALOpen()
  in order to do this, but I'm not exactly sure how.  Could someone pass me
  a
  couple lines of Python that show how this is done?
 
  Thanks in advance.
 
  Roger
  --
  ___
  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] gdalwarp interpolation artifacts

[Roger André]

Hi Seth and Frank,

Thanks for the feedback.  I just finished reading the changelog link that
Seth sent and all I can say is, Wow, nice job!

I did some further testing last night using a single GTOPO30 tile, and doing
a 10x upsample.  The artifacts are present in it at columns 6000 and 12000,
so I think we can rule out that gdal_merge.py is introducing the error.  I'm
testing the cubic resampling operator this morning, in the hopes that it
won't produce the error.  Since it seems to require a very large scaling
difference, it also takes a bit of time to run the test.  I'll let you know
what I find when it's done.

Regarding running the latest code from trunk, I'm not sure if I can do
that.  I'm currently using the gdal implementation in the Linux version of
FWTools-2.0.6 because I am unable to get all of the dependencies installed
that I need for a source-compiled version of gdal to work properly on my
system.  I think I remember seeing some traffic from Mateuz in the recent
past that indicated he used FWTools to build against.  I'll see if I can
figure out what he was doing.  If I recall correctly, Frank wasn't a huge
fan of this technique though. ;)

Thanks again,

Roger
--

On Mon, Sep 29, 2008 at 10:34 PM, Seth Price [EMAIL PROTECTED] wrote:

 I agree with Frank, this sounds like the bug that I reported (and fixed)
 here:
 http://trac.osgeo.org/gdal/ticket/2327

 (Although it became a log of my alterations to the warper kernel, which was
 well beyond the scope of the original bug.)
 ~Seth



 On Sep 29, 2008, at 11:19 PM, Frank Warmerdam wrote:

  On Tue, Sep 30, 2008 at 4:30 AM, Roger André [EMAIL PROTECTED] wrote:

 Hi List,

 I've hit the same problem twice now, and I'm pretty certain after Round 2
 that I'm not introducing it.  Basically, I'm using gdal_merge.py to
 mosaic a
 group of low-resolution, 32-bit floating point rasters together, then
 running gdalwarp -ts xxx yyy -r cubicspline on the mosaic to get a
 higher
 resolution image.  I'm finding afterwards that there are linear artifacts
 in
 the mosaic which run vertically through the mosaic, which look like edge
 artifacts, but which are not located near the edge of either an original
 low-res tile, or the resulting high res one.  I've tested the
 interpolation
 of a single low-res tile in the area where the artifact is present in the
 mosaic and I do not get the same results - the interpolated tile is clear
 of
 artifacts.

 Is it possible that I'm hitting some sort of memory limit that is causing
 only a certain number of columns to be interpolated at one time, rather
 than
 the entire row?


 Roger,

 There have been some fixes to some of the gdalwarp interpolator
 code in recent months, so one hint is to ensure you try the latest
 trunk code to see if perhaps the problem is already fixed.

 Second, you may want to investigate the gdal_merge.py product
 carefully to ensure there aren't any bad pixels in it along the
 boundaries.

 Third you might try the cubic (instead of cubicspline) interpolator
 to see if it gives cleaner results.

 Beyond that a ticket submitted on the smallest test case that
 demonstrates the problem would be appreciated.

 Best regards,
 --

 ---+--
 I set the clouds in motion - turn up   | Frank Warmerdam,
 [EMAIL PROTECTED]
 light and sound - activate the windows | 
 http://pobox.com/~warmerdamhttp://pobox.com/%7Ewarmerdam
 and watch the world go round - Rush| Geospatial Programmer for Rent
 ___
 gdal-dev mailing list
 gdal-dev@lists.osgeo.org
 http://lists.osgeo.org/mailman/listinfo/gdal-dev



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

Re: [gdal-dev] creating an epsg code that can be used within the epsg file for fwtools and MapServer

[Roger André]

Hi John,

No guarantees that I got this right, but you could try adding the following
entry to your FWTools/share/proj/epsg file:

# Guatemala Transverse Mercator
001234 +proj=tmerc +x_0=50.00 +y_0=0 +lon_0=-90.5 +lat_0=0
+k=0.999800 +ellps=WGS84 +datum=WGS84 +to_meter=0.3048006096012192 +no_defs


Then you should be able to acces it by using EPSG:001234 in your -s_srs
parameters.

Roger
--

On Mon, Oct 6, 2008 at 7:10 AM, John Mitchell [EMAIL PROTECTED]wrote:

 Hi,

 How would I create an epsg code that can be used within the epsg file for
 fwtools and MapServer?

 Below is the information that I would like to use to create an epsg code.
 What would be format for the new epsg code?  Is the information listed below
 enough to create an epsg code and if not what am I missing?


 *The 2001 data was originally in a Guatemalan Transverse Mercator
 coordinate
 system with the parameters shown below:

 Projection: Transverse Mercator
 False Easting: 50.00
 False Northing: 0.00
 Central Meridian: -90.50
 Scale Factor: 0.999800
 Latitude of Origin: 0.000
 Linear Unit: Meter
 Geographic Coordinate System: GCS_WGS_1984

 *Below is the gdalinfo listing for one of files that I would like to
 associate to this new epsg code:

 *Driver: GTiff/GeoTIFF
 Size is 11214, 7860
 Coordinate System is `'
 Origin = (317049.000,**1659460.000)
 Pixel Size = (0.500,-0.**500)
 Metadata:
   AREA_OR_POINT=Area
   TIFFTAG_SOFTWARE=ERDAS IMAGINE
   TIFFTAG_XRESOLUTION=1
   TIFFTAG_YRESOLUTION=1
   TIFFTAG_RESOLUTIONUNIT=1 (unitless)
 Corner Coordinates:
 Upper Left  (  317049.000, 1659460.000)
 Lower Left  (  317049.000, 1655530.000)
 Upper Right (  322656.000, 1659460.000)
 Lower Right (  322656.000, 1655530.000)
 Center  (  319852.500, 1657495.000)
 Band 1 Block=11214x8 Type=Byte, ColorInterp=Red
 Band 2 Block=11214x8 Type=Byte, ColorInterp=Green
 Band 3 Block=11214x8 Type=Byte, ColorInterp=Blue*



 Thanks,

 --
 John J. Mitchell

 ___
 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] nodata value assigned via script

[Roger André]

Hi,

I'm getting raster data written by different people who can't seem to use
the same nodata values.  I'm reading the data files via the Python API to
GDAL, and doing various calculations on the data values, then writing new
rasters with the results of my calculations.  Is there a nodata value that
I can assign to a pixel that the gdal utilities will recognize as such later
on?  All my other values are 32-bit floating point values, so I'm not sure
what value to assign.

Thanks.

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

[gdal-dev] Problem reprojecting from EPSG:4326 to EPSG:900913

[Roger André]

Hi Everyone,

I'm having some issues when reprojecting data via a WMS request to
MapServer.  While I realize this only touches on the gdal side of things
slightly (I'm using the gdal driver and FWTools version 2.0.6 in Linux), I'm
hoping the problem is one that someone on this list may have encountered.
Very briefly, here's what's happening.  I have an image file that is
completely unprojected, and which spans a global extent in pixels that are
.025 x .025 degrees in size.  I have assigned this file an SRS of EPSG:4326,
and am trying to display it with an SRS of 900913 in Google Maps.  The
reprojection seems to work, at least in that the edges of the data set line
up correctly with all of the features in Google Maps, but there are
horizontal lines across my data in various places.  When I replicate the WMS
call for one of these bad tiles,

http://localhost/cgi-bin/mapserv?
map=/var/www/mapfiles/test_20m/test.map
layers=imagery
srs=EPSG%3A900913
version=1.1.1
bbox=-20037508.34%2C-20037508.34%2C-10018754.1704%2C-10018754.1704
service=WMS
width=256
styles=
format=image%2Fpng
request=GetMap
height=256

I get back an image that is 256 x 256 pixels, but where the bottom 4 rows of
pixels are essentially empty (they're filled with the MapServer background
color).  The pixels that are correctly represented are also stretched a bit
in the Y-axis as well.  If I screw around with my projection definition for
EPSG:900913 and add an ellipsoid definition, +ellps=WGS84, the missing 4
pixels will be correctly rendered in this particular example, but if I look
closely at other parts of the data, in other zoom levels, the same problem
reappears.

I've also tried reprojecting the data to EPSG:900913 first, and then
displaying the data in Google Maps without making any sort of reprojection
request.  What I get then from the image request above is a tile that is
filled 2/3 with correctly sized pixels (no y-axis stretch), and 1/3 is
empty.  So I'm guessing that my definition of Google's projection is not
correct, but I have no clue how to fix it.  I've used the definition found
at http://spatialreference.org/ref/user/google-projection/, and also at
http://trac.openlayers.org/wiki/SphericalMercator.

900913 +proj=merc +ellps=WGS84 +a=6378137 +b=6378137 +lat_ts=0.0
+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m [EMAIL PROTECTED] +no_defs 

It seems pretty well researched at both of these sites, but I have never
seen the [EMAIL PROTECTED] parameter before.  I've read the info about it at
http://proj.maptools.org/faq.html, but don't fully understand it.

Thanks in advance for any help or advice you can give me.

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

[gdal-dev] Going between gdal nodata and GMT NaN

[Roger André]

Sooo... apparently nodata values are not converted to NaN when doing
something like a gdal_translate from GeoTiff to GMT format.  Is there a way
that anyone knows to do this?

Thanks,

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

[gdal-dev] Re: Going between gdal nodata and GMT NaN

[Roger André]

Hmm, let me ask this a different way.  Does the GMT driver implemented in
gmtdataset.cpp support the creation of Nan values when writing to an
output dataset?

Thanks,

Roger
--

On Wed, Oct 22, 2008 at 4:38 PM, Roger André [EMAIL PROTECTED] wrote:

 Sooo... apparently nodata values are not converted to NaN when doing
 something like a gdal_translate from GeoTiff to GMT format.  Is there a way
 that anyone knows to do this?

 Thanks,

 Roger
 --

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

Re: [gdal-dev] gdal_rasterize usage problems

[Roger André]

Hi Markus,

A friend of mine just tried this out on r15709 and reports that it works
fine.  Are you referring to the syntax in the original email not working, or
what I advised Andrew to try?  If it's the latter, please send the command
you're using, so we can dissect it.

Thanks,

Roger
--

On Tue, Nov 11, 2008 at 2:08 AM, Markus Neteler [EMAIL PROTECTED] wrote:

 On Mon, Nov 10, 2008 at 10:43 PM, Roger André [EMAIL PROTECTED] wrote:
  Andrew,
 
  I believe your syntax is incorrect.  Assuming that you are trying to
 exclude
  water from your image, it should be something like this:
  % gdal_rasterize -b 1 -b 2 -b 3 -i -burn 0 -burn 0 -burn 0 -l gshhs_land
  gshhs_land.shp warped.tif
 
  Have a look at http://www.gdal.org/gdal_rasterize.html for a list of all
 the
  options.

 Roger,

 it fails also for me (using recent trunk).

 Markus

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

[gdal-dev] Can I change SOURCE_EXTRA setting at compile time?

[Roger André]

I've been mercilessly beating a dead horse on a couple of lists about some
weird artifacts that I'm getting when doing reprojections in Mapserver from
epsg:4326 to epsg:900913.  It's recently come to my attention that something
was discussed in an IRC chat that might be pertinent to this issue, and I
was hoping I could test it.  Specifically, it relates to SOURCE_EXTRA, and
although the conversation at the time dealt with issues at the International
dateline, I think it might also manifest itself in other edge cases (literal
edges). ;).

I would like to know if it is possible to set SOURCE_EXTRA to something
other than 1 at compile time?  I do not see it listed as an option in the
configure documentation, but perhaps I could edit a file, or pass in the
flag anyhow?

Thanks in advance,

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

Re: [gdal-dev] Can I change SOURCE_EXTRA setting at compile time?

[Roger André]

Hi Frank,

I've theorized all sorts of reasons for what I'm seeing.  So far none of
them are right.  Basically, if I take a really high-res epsg:4326 image, and
serve it up through Mapserver in epsg:900913, there are instances where I
get images back which are missing rows of pixels at the bottom of the
image.  If I reduce the resolution of the image the lines go away.  I
haven't yet figured out at what resolution the problem starts to appear.

I just tested the idea of changing the values in the warper code:
First:
nResWinSize += 100;
Then:
nResWinSize += 500;

Then rebuilt gdal.   Didn't seem to make a difference.

Thanks for the info though.

Roger
--

On Tue, Nov 18, 2008 at 5:07 PM, Frank Warmerdam [EMAIL PROTECTED]wrote:

 Roger André wrote:

 I've been mercilessly beating a dead horse on a couple of lists about some
 weird artifacts that I'm getting when doing reprojections in Mapserver from
 epsg:4326 to epsg:900913.  It's recently come to my attention that something
 was discussed in an IRC chat that might be pertinent to this issue, and I
 was hoping I could test it.  Specifically, it relates to SOURCE_EXTRA, and
 although the conversation at the time dealt with issues at the International
 dateline, I think it might also manifest itself in other edge cases (literal
 edges). ;).

 I would like to know if it is possible to set SOURCE_EXTRA to something
 other than 1 at compile time?  I do not see it listed as an option in the
 configure documentation, but perhaps I could edit a file, or pass in the
 flag anyhow?


 Roger,

 I haven't followed your problems closely (I recall something about datum
 shifts).  The SOURCE_EXTRA only helps with the GDAL warper.  If you wanted
 you could modify the following code in gdal/alg/gdalwarpoperation.cpp
 and recompile:

 /*  */
 /*  Allow addition of extra sample pixels to source window to   */
 /*  avoid missing pixels due to sampling error.  In fact,   */
 /*  fallback to adding a bit to the window if any points failed */
 /*  to transform.   */
 /*  */
if( CSLFetchNameValue( psOptions-papszWarpOptions,
   SOURCE_EXTRA ) != NULL )
{
nResWinSize += atoi(
CSLFetchNameValue( psOptions-papszWarpOptions, SOURCE_EXTRA
 ));
}
else if( nFailedCount  0 )
nResWinSize += 10;


 Perhaps change the nResWinSize += 10 line to nResWinSize += 100 or
 something
 like that.

 Best regards,
 --

 ---+--
 I set the clouds in motion - turn up   | Frank Warmerdam,
 [EMAIL PROTECTED]
 light and sound - activate the windows | 
 http://pobox.com/~warmerdamhttp://pobox.com/%7Ewarmerdam
 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] Can I change SOURCE_EXTRA setting at compile time?

[Roger André]

Hi Frank,

You're absolutely right.  Changing those values in mapresample.c to 10
removed the artifacts when using a WMS request.  They're still present when
the request is made via mapscript though.  I need to manually code a
mapscript example though, because right now I'm testing that via TileCache,
and I'm not 100% certain how the image is being requested.  Can you think of
a reason why WMS would be different from a straight read and reproject
through mapscript?

cc'ing TileCache list and mapserver lists as well
--

On Tue, Nov 18, 2008 at 5:58 PM, Frank Warmerdam [EMAIL PROTECTED]wrote:

 Roger André wrote:

 Hi Frank,

 I've theorized all sorts of reasons for what I'm seeing.  So far none of
 them are right.  Basically, if I take a really high-res epsg:4326 image, and
 serve it up through Mapserver in epsg:900913, there are instances where I
 get images back which are missing rows of pixels at the bottom of the image.
  If I reduce the resolution of the image the lines go away.  I haven't yet
 figured out at what resolution the problem starts to appear.

 I just tested the idea of changing the values in the warper code:
 First:
 nResWinSize += 100;
 Then:
 nResWinSize += 500;

 Then rebuilt gdal.   Didn't seem to make a difference.


 Roger,

 Are you sure the GDAL warper is being used? If the reprojection is being
 done by MapServer it will be using the MapServer warper.

 Try changing the +- 1.0 values to larger values in the following code in
 mapserver's mapresample.c:

 /*  */
 /*  Project desired extents out by 2 pixels, and then strip to  */
 /*  available data. */
 /*  */
memcpy( sOrigSrcExtent, sSrcExtent, sizeof(sSrcExtent) );

sSrcExtent.minx = floor(sSrcExtent.minx-1.0);
sSrcExtent.maxx = ceil (sSrcExtent.maxx+1.0);
sSrcExtent.miny = floor(sSrcExtent.miny-1.0);
sSrcExtent.maxy = ceil (sSrcExtent.maxy+1.0);


 See if that helps.  Currently there is no processing option that does
 what SOURCE_EXTRA does, though it might be worth adding one.


 Best regards,
 --

 ---+--
 I set the clouds in motion - turn up   | Frank Warmerdam,
 [EMAIL PROTECTED]
 light and sound - activate the windows | 
 http://pobox.com/~warmerdamhttp://pobox.com/%7Ewarmerdam
 and watch the world go round - Rush| Geospatial Programmer for Rent


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

[gdal-dev] Re: gdal segmentation fault when reading NetCDF grid files

[Roger André]

Sorry, should have trolled the Web first.  I see now that this is nothing
new, and is caused by a conflict between libnetcdf and libhdf4/5.  I've
resolved this by doing the following:

1. Removing all hdf libs from my system:
$ sudo apt-get remove libhdf4g libhdf4g-dev libhdf5-serial-1.6.5-0
libhdf5-serial-dev

2. Doing a sudo make clean in my /usr/local/gdal dir

3. Doing a new configure with --with-hdf4=no --with-hdf5=no, then running
make.

Seems to work now,

Thanks.
--

On Mon, Nov 24, 2008 at 10:07 AM, Roger André [EMAIL PROTECTED] wrote:

 Hi everyone,

 I'm having some problems getting gdal to read netCDF grid files using
 gdal-1.5.2 which I have built from source.  It is running on an Ubuntu 8.04
 machine, and was built after the following dependencies were installed:

 $ sudo apt-get install libnetcdf4 libnetcdf-dev netcdfg-dev netcdf-doc
 libhdf4g libhdf4g-dev libhdf5-serial-1.6.5-0 libhdf5-serial-1.6.5-0
 libhdf5-serial-dev

 I built gdal with the following options:

 $ sudo ./configure --with-libz --with-pg=/usr/local/pgsql/bin/pg_config
 --with-netcdf --with-png --with-geotiff --with-jpeg --with-gif --with-hdf4
 --with-hdf5 --with-jasper --with-geos=/usr/local/bin/geos-config
 --with-python

 When I try to open the file with gdalinfo, I get the following seg fault in
 /var/log/syslog:

 Nov 24 09:42:34 randre-linux kernel: [ 1921.499150] gdalinfo[6971]:
 segfault at 1031 eip b713d886 esp bfefb390 error 4

 The same file can be opened with no problem using gdalinfo from
 FWTools-2.0.6

 Is it possible I am missing some netCDF/HDF fu needed for gdal to work with
 these sorts of files?

 Thanks in advance for any suggestions,

 Roger
 --



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

[gdal-dev] Working between GeoTIFF's and Arc/Info ASCII Grid formats with NAN values

[Roger André]

Hi All,

Just had an email from a client complaining about the nan values that I
included in some Arc/Info ASCII Grid files that I sent them.  The files were
created by converting GeoTIFF's into ASCII Grid format via gdal_translate.
I personally like that the nan's were carried through, but I would like to
know if there is a way that I can map the nan into a - value via the
gdal utilities?

Thanks,

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

Re: [gdal-dev] Working between GeoTIFF's and Arc/Info ASCII Grid formats with NAN values

[Roger André]

ROTFL  You could also potentially do it with by ascii search and replace in
the resulting grid file

You are SO right!  I had totally forgotten that this was ASCII output.  Doh!

Thanks Frank.
--



On Tue, Dec 9, 2008 at 5:24 PM, Frank Warmerdam [EMAIL PROTECTED] wrote:

 Roger André wrote:

 Hi All,

 Just had an email from a client complaining about the nan values that I
 included in some Arc/Info ASCII Grid files that I sent them.  The files were
 created by converting GeoTIFF's into ASCII Grid format via gdal_translate.
  I personally like that the nan's were carried through, but I would like
 to know if there is a way that I can map the nan into a - value via
 the gdal utilities?


 Roger,

 I'm not aware of a way to convert NaN values to an actual value with the
 GDAL command line utilities.  You could accomplish this with a program,
 or a python script.  You could also potentially do it with by ascii
 search and replace in the resulting grid file.

 Best regards,
 --

 ---+--
 I set the clouds in motion - turn up   | Frank Warmerdam,
 [EMAIL PROTECTED]
 light and sound - activate the windows | 
 http://pobox.com/~warmerdamhttp://pobox.com/%7Ewarmerdam
 and watch the world go round - Rush| Geospatial Programmer for Rent


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

[gdal-dev] gdal2tiles.py and dateline gap

[Roger André]

I've gotten the map tile output from gdal2tiles.py to work properly with
Tilecache and the results look pretty good.  However, I'm seeing a fairly
large gap in data at the international dateline.  I recall seeing some
discussion about this between Frank and Klokan on the list previously, and
recompiled gdal after editing gdal/alg/gdalwarpoperation.cpp and adding a
larger size to the sampling edge buffer.  This doesn't seem to have helped
though.

Is it possible to create a global tileset with gdal2tiles.py without having
a gap at the dateline?

Thanks,

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

Re: [gdal-dev] gdal_translate input separate R,G,B output combined RGB

[Roger André]

Hi Greg,

I don't know of a gdal_translate command for this, but you can use the
gdal_merge.py tool to do this.

gdal_merge.py - o outfile -separate gtif_R.tif gtif_G.tif gtif_B.tif

Roger
--

2009/1/9 Greg Coats gregco...@mac.com

 I have 3 single band GeoTIF images of the same geographic area: one for the
 Red band, one for the Green band, and one for the Blue
 band, gtif_R.tif, gtif_G.tif, gtif_B.tif.  Please, what is the
 gdal_translate syntax http://www.gdal.org/gdal_translate.html 
 http://www.gdal.org/gdal_translate.htmlto
 combine these 3 single band GeoTIF images into 1 three band RGB GeoTIF
 image, gtif_RGB.tif?  Greg

 ___
 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] Trying to clean up multiple GDAL installs

[Roger André]

Ok, I admit it.  I was lured by the promise of what I didn't already have.
I had a perfectly functional gdal install from source on my Ubuntu box.
Then I went and installed Qgis using apt-get, which brought in a gdal
package as well, and started to get errors like the one below.

ogrinfo: symbol lookup error: ogrinfo: undefined symbol:
_ZN10OGRFeature12DumpReadableEP8_IO_FILE

I decided I'd rather have ogrinfo than Qgis, so I removed Qgis using
apt-get autoremove qgis.  This said that it removed some gdal stuff, but
afterwards I still had the ogrinfo error.  I then went into my
/usr/local/gdal-1.5.2/ and did a make clean.  I also deleted all ogr
and gdal apps that were left in /usr/local/bin.  I then built gdal from
source again in the /usr/local/gdal-1.5.2/ dir.  However, when I do a
ogrinfo --version it tells me I'm using GDAL 1.6.0 - and I still have
the ogrinfo error.

I've checked my gdal-config file and it says CONFIG_VERSION=1.5.2, which
is correct.  So now I'm kinda scratching my head wondering what to delete
next.

Any hints, please?

Thanks,

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

[gdal-dev] NetCDF data upside down when read by gdal

[Roger André]

Hi All,

- I am attempting to use a NetCDF file as a data source for a raster layer
in Mapserver.  According to a very useful email sent in November 2005 from
Norman Barker,
http://lists.osgeo.org/pipermail/mapserver-users/2005-November/012153.html,
and from Mapserver documentation at
http://www.mapserver.org/ogc/wcs_format.html#netcdf, it seems like this
should be possible using Mapserver's gdal driver to read the file.  However,
I am running into some problems when trying to implement these examples
using GDAL 1.5.2, released 2008/05/29.

- The NetCDF file which I am using to test with is a CF Convention sample
file which I downloaded from
http://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc.
This file contains information about sea-surface temperatures, and when
viewed by gdalinfo has the following characteristics:

$ gdalinfo tos_O1_2001-2002.nc
Driver: netCDF/Network Common Data Format
Files: tos_O1_2001-2002.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#title=IPSL  model output prepared for IPCC Fourth Assessment
SRES A2 experiment
  NC_GLOBAL#institution=IPSL (Institut Pierre Simon Laplace, Paris, France)
  NC_GLOBAL#source=IPSL-CM4_v1 (2003) : atmosphere : LMDZ (IPSL-CM4_IPCC,
96x71x19) ; ocean ORCA2 (ipsl_cm4_v1_8, 2x2L31); sea ice LIM (ipsl_cm4_v
  NC_GLOBAL#contact=Sebastien Denvil, sebastien.den...@ipsl.jussieu.fr
  NC_GLOBAL#project_id=IPCC Fourth Assessment
  NC_GLOBAL#table_id=Table O1 (13 November 2004)
  NC_GLOBAL#experiment_id=SRES A2 experiment
  NC_GLOBAL#realization=1
  NC_GLOBAL#cmor_version=9.60e-01
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#history=/MM/JJ: data generated; /MM/JJ+1 data
transformed  At 16:37:23 on 01/11/2005, CMOR rewrote data to comply with CF
standards and IPCC Fourth Assessment requirements
  NC_GLOBAL#references=Dufresne et al, Journal of Climate, 2015, vol XX, p
136
  NC_GLOBAL#comment=Test drive
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:tos_O1_2001-2002.nc:lon_bnds
  SUBDATASET_1_DESC=[180x2] lon_bnds (64-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:tos_O1_2001-2002.nc:lat_bnds
  SUBDATASET_2_DESC=[170x2] lat_bnds (64-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:tos_O1_2001-2002.nc:time_bnds
  SUBDATASET_3_DESC=[24x2] time_bnds (64-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:tos_O1_2001-2002.nc:tos
  SUBDATASET_4_DESC=[24x170x180] sea_surface_temperature (32-bit
floating-point)
Corner Coordinates:
Upper Left  (0.0,0.0)
Lower Left  (0.0,  512.0)
Upper Right (  512.0,0.0)
Lower Right (  512.0,  512.0)
Center  (  256.0,  256.0)

- The SUBDATASET of interest is the tos one, and when queried by gdalinfo
it shows that it contains 24 bands, each of which varies by
NETCDF_DIMENSION_time that increments by 15.  When specifically queried by
gdalinfo, this SUBDATASET looks like this:

$ gdalinfo NETCDF:tos_O1_2001-2002.nc:tos
Driver: netCDF/Network Common Data Format
Files: none associated
Size is 180, 170
Coordinate System is `'
Origin = (0.000,90.000)
Pixel Size = (2.000,-1.000)
Metadata:
  NC_GLOBAL#title=IPSL  model output prepared for IPCC Fourth Assessment
SRES A2 experiment
  NC_GLOBAL#institution=IPSL (Institut Pierre Simon Laplace, Paris, France)
  NC_GLOBAL#source=IPSL-CM4_v1 (2003) : atmosphere : LMDZ (IPSL-CM4_IPCC,
96x71x19) ; ocean ORCA2 (ipsl_cm4_v1_8, 2x2L31); sea ice LIM (ipsl_cm4_v
  NC_GLOBAL#contact=Sebastien Denvil, sebastien.den...@ipsl.jussieu.fr
  NC_GLOBAL#project_id=IPCC Fourth Assessment
  NC_GLOBAL#table_id=Table O1 (13 November 2004)
  NC_GLOBAL#experiment_id=SRES A2 experiment
  NC_GLOBAL#realization=1
  NC_GLOBAL#cmor_version=9.60e-01
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#history=/MM/JJ: data generated; /MM/JJ+1 data
transformed  At 16:37:23 on 01/11/2005, CMOR rewrote data to comply with CF
standards and IPCC Fourth Assessment requirements
  NC_GLOBAL#references=Dufresne et al, Journal of Climate, 2015, vol XX, p
136
  NC_GLOBAL#comment=Test drive
  tos#standard_name=sea_surface_temperature
  tos#long_name=Sea Surface Temperature
  tos#units=K
  tos#cell_methods=time: mean (interval: 30 minutes)
  tos#_FillValue=1.00e+20
  tos#missing_value=1.00e+20
  tos#original_name=sosstsst
  tos#original_units=degC
  tos#history= At   16:37:23 on 01/11/2005: CMOR altered the data in the
following ways: added 2.73150E+02 to yield output units;  Cyclical dimension
was output starting at a different lon;
  lon#standard_name=longitude
  lon#long_name=longitude
  lon#units=degrees_east
  lon#axis=X
  lon#bounds=lon_bnds
  lon#original_units=degrees_east
  lat#standard_name=latitude
  lat#long_name=latitude
  lat#units=degrees_north
  lat#axis=Y
  lat#bounds=lat_bnds
  lat#original_units=degrees_north
  time#standard_name=time
  time#long_name=time
  time#units=days since 2001-1-1
  time#axis=T
  time#calendar=360_day
  time#bounds=time_bnds
  time#original_units=seconds since 2001-1-1
Corner Coordinates:

Re: [gdal-dev] NetCDF data upside down when read by gdal

[Roger André]

Rockin`!

Thanks so much, that helps very much.

Roger
--

On Fri, Feb 20, 2009 at 8:32 AM, Gong, Shawn (Contractor) 
shawn.g...@drdc-rddc.gc.ca wrote:

  hi Roger,



 Here are cpp and h files that I modified to flip netCDF images.

 They are in ..\gdal-1.6.0\frmts\netcdf\



 hope that they help,
 Shawn


  --

 *From:* gdal-dev-boun...@lists.osgeo.org [mailto:
 gdal-dev-boun...@lists.osgeo.org] *On Behalf Of *Roger André
 *Sent:* Wednesday, February 18, 2009 3:45 PM
 *To:* gdal-dev
 *Subject:* [gdal-dev] NetCDF data upside down when read by gdal



 Hi All,

 - I am attempting to use a NetCDF file as a data source for a raster layer
 in Mapserver.  According to a very useful email sent in November 2005 from
 Norman Barker,
 http://lists.osgeo.org/pipermail/mapserver-users/2005-November/012153.html,
 and from Mapserver documentation at
 http://www.mapserver.org/ogc/wcs_format.html#netcdf, it seems like this
 should be possible using Mapserver's gdal driver to read the file.  However,
 I am running into some problems when trying to implement these examples
 using GDAL 1.5.2, released 2008/05/29.

 - The NetCDF file which I am using to test with is a CF Convention sample
 file which I downloaded from
 http://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc.
 This file contains information about sea-surface temperatures, and when
 viewed by gdalinfo has the following characteristics:

 $ gdalinfo tos_O1_2001-2002.nc
 Driver: netCDF/Network Common Data Format
 Files: tos_O1_2001-2002.nc
 Size is 512, 512
 Coordinate System is `'
 Metadata:
   NC_GLOBAL#title=IPSL  model output prepared for IPCC Fourth Assessment
 SRES A2 experiment
   NC_GLOBAL#institution=IPSL (Institut Pierre Simon Laplace, Paris, France)
   NC_GLOBAL#source=IPSL-CM4_v1 (2003) : atmosphere : LMDZ (IPSL-CM4_IPCC,
 96x71x19) ; ocean ORCA2 (ipsl_cm4_v1_8, 2x2L31); sea ice LIM (ipsl_cm4_v
   NC_GLOBAL#contact=Sebastien Denvil, sebastien.den...@ipsl.jussieu.fr
   NC_GLOBAL#project_id=IPCC Fourth Assessment
   NC_GLOBAL#table_id=Table O1 (13 November 2004)
   NC_GLOBAL#experiment_id=SRES A2 experiment
   NC_GLOBAL#realization=1
   NC_GLOBAL#cmor_version=9.60e-01
   NC_GLOBAL#Conventions=CF-1.0
   NC_GLOBAL#history=/MM/JJ: data generated; /MM/JJ+1 data
 transformed  At 16:37:23 on 01/11/2005, CMOR rewrote data to comply with CF
 standards and IPCC Fourth Assessment requirements
   NC_GLOBAL#references=Dufresne et al, Journal of Climate, 2015, vol XX, p
 136
   NC_GLOBAL#comment=Test drive
 Subdatasets:
   SUBDATASET_1_NAME=NETCDF:tos_O1_2001-2002.nc:lon_bnds
   SUBDATASET_1_DESC=[180x2] lon_bnds (64-bit floating-point)
   SUBDATASET_2_NAME=NETCDF:tos_O1_2001-2002.nc:lat_bnds
   SUBDATASET_2_DESC=[170x2] lat_bnds (64-bit floating-point)
   SUBDATASET_3_NAME=NETCDF:tos_O1_2001-2002.nc:time_bnds
   SUBDATASET_3_DESC=[24x2] time_bnds (64-bit floating-point)
   SUBDATASET_4_NAME=NETCDF:tos_O1_2001-2002.nc:tos
   SUBDATASET_4_DESC=[24x170x180] sea_surface_temperature (32-bit
 floating-point)
 Corner Coordinates:
 Upper Left  (0.0,0.0)
 Lower Left  (0.0,  512.0)
 Upper Right (  512.0,0.0)
 Lower Right (  512.0,  512.0)
 Center  (  256.0,  256.0)

 - The SUBDATASET of interest is the tos one, and when queried by gdalinfo
 it shows that it contains 24 bands, each of which varies by
 NETCDF_DIMENSION_time that increments by 15.  When specifically queried by
 gdalinfo, this SUBDATASET looks like this:

 $ gdalinfo NETCDF:tos_O1_2001-2002.nc:tos
 Driver: netCDF/Network Common Data Format
 Files: none associated
 Size is 180, 170
 Coordinate System is `'
 Origin = (0.000,90.000)
 Pixel Size = (2.000,-1.000)
 Metadata:
   NC_GLOBAL#title=IPSL  model output prepared for IPCC Fourth Assessment
 SRES A2 experiment
   NC_GLOBAL#institution=IPSL (Institut Pierre Simon Laplace, Paris, France)
   NC_GLOBAL#source=IPSL-CM4_v1 (2003) : atmosphere : LMDZ (IPSL-CM4_IPCC,
 96x71x19) ; ocean ORCA2 (ipsl_cm4_v1_8, 2x2L31); sea ice LIM (ipsl_cm4_v
   NC_GLOBAL#contact=Sebastien Denvil, sebastien.den...@ipsl.jussieu.fr
   NC_GLOBAL#project_id=IPCC Fourth Assessment
   NC_GLOBAL#table_id=Table O1 (13 November 2004)
   NC_GLOBAL#experiment_id=SRES A2 experiment
   NC_GLOBAL#realization=1
   NC_GLOBAL#cmor_version=9.60e-01
   NC_GLOBAL#Conventions=CF-1.0
   NC_GLOBAL#history=/MM/JJ: data generated; /MM/JJ+1 data
 transformed  At 16:37:23 on 01/11/2005, CMOR rewrote data to comply with CF
 standards and IPCC Fourth Assessment requirements
   NC_GLOBAL#references=Dufresne et al, Journal of Climate, 2015, vol XX, p
 136
   NC_GLOBAL#comment=Test drive
   tos#standard_name=sea_surface_temperature
   tos#long_name=Sea Surface Temperature
   tos#units=K
   tos#cell_methods=time: mean (interval: 30 minutes)
   tos#_FillValue=1.00e+20
   tos#missing_value=1.00e+20
   tos#original_name=sosstsst
   tos#original_units=degC
   tos#history= At   16

[gdal-dev] gdal, netCDF and projection

[Roger André]

Hi all,

Just a couple quick netCDF related questions:

1)  Is there any way to include projection information in a netCDF
file so that gdal can properly read it?  It appears that
georeferencing is handled correctly, but not projection info.  This
causes problems with the whole north/south origin shift between the
netCDF and other raster formats when gdal is used to convert the .nc
file into a GeoTIFF.  Most notably, I'm not sure how to assign the
data a projection - either for display in MapServer, or to reproject
it to EPSG:4326:  See example below.

File created with ArcGIS 9.3:
$ gdalinfo fwdir_bas.nc
Driver: netCDF/Network Common Data Format
Files: fwdir_bas.nc
Size is 42000, 42000
Coordinate System is `'
Origin = (65.000,5.0028436)
Pixel Size = (0.0008333,0.0008333)
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#Source_Software=ESRI ArcGIS
  Fwdir_bas#long_name=Fwdir_bas
 
Fwdir_bas#esri_pe_string=GEOGCS[GCS_WGS_1984,DATUM[D_WGS_1984,SPHEROID[WGS_1984,6378137.0,298.257223563]],PRIMEM[Greenwich,0.0],UNIT[Degree,0.0174532925199433]]
  Fwdir_bas#coordinates=lon lat
  Fwdir_bas#units=Degree
  Fwdir_bas#missing_value=
  lon#long_name=longitude coordinate
  lon#standard_name=longitude
  lon#units=degrees_east
  lat#long_name=latitude coordinate
  lat#standard_name=latitude
  lat#units=degrees_north
Corner Coordinates:
Upper Left  (  65.000,   5.000)
Lower Left  (  65.000,  40.000)
Upper Right ( 100.000,   5.000)
Lower Right ( 100.000,  40.000)
Center  (  82.500,  22.500)
Band 1 Block=42000x1 Type=Byte, ColorInterp=Undefined
  NoData Value=0
  Metadata:
NETCDF_VARNAME=Fwdir_bas

- After translation to GeoTIFF:

$ gdalinfo fwdir_bas.tif
Driver: GTiff/GeoTIFF
Files: fwdir_bas.tif
Size is 42000, 42000
Coordinate System is `'
Origin = (65.000,5.0028436)
Pixel Size = (0.0008333,0.0008333)
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#Source_Software=ESRI ArcGIS
  Fwdir_bas#long_name=Fwdir_bas
  
Fwdir_bas#esri_pe_string=GEOGCS[GCS_WGS_1984,DATUM[D_WGS_1984,SPHEROID[WGS_1984,6378137.0,298.257223563]],PRIMEM[Greenwich,0.0],UNIT[Degree,0.0174532925199433]]
  Fwdir_bas#coordinates=lon lat
  Fwdir_bas#units=Degree
  lon#long_name=longitude coordinate
  lon#standard_name=longitude
  lon#units=degrees_east
  lat#long_name=latitude coordinate
  lat#standard_name=latitude
  lat#units=degrees_north
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  65.000,   5.000)
Lower Left  (  65.000,  40.000)
Upper Right ( 100.000,   5.000)
Lower Right ( 100.000,  40.000)
Center  (  82.500,  22.500)
Band 1 Block=42000x1 Type=Byte, ColorInterp=Gray
  NoData Value=0
  Metadata:
NETCDF_VARNAME=Fwdir_bas

2)  Is there any continuing work being done in the gdal community with
regards to the GALEON project?  As far as I can see, all gdal posts
related to this are from 2005 - 2006, yet recent articles in some
industry rags suggest that GALEON is still on the OGC's radar.  I'd be
interested in working on this with anyone else who already is to help
resolve some of the remaining hiccups that are left in getting netCDF
files to work more easily with gdal.  Feel free to contact me off-list
regarding this.

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

Re: [gdal-dev] gdal, netCDF and projection

[Roger André]

Hi Frank,

Sorry for being unclear about the projection stuff.  I guess what I
was asking was how to go about defining the South-up nature of the
data.  But I guess it doesn't matter, so long as the corner coords and
pixel dimensions are correct.  I used gdal_translate -a_srs
EPSG:4326 on the GeoTIFF I created earlier, and it seems to have
assigned the SRS properly and kept a positive Y axis pixel size.  This
South-Up stuff is just weird to me, I guess.

$ gdalinfo wgs84_fwdir_bas.tif
Driver: GTiff/GeoTIFF
Files: wgs84_fwdir_bas.tif
Size is 42000, 42000
Coordinate System is:
GEOGCS[WGS 84,
DATUM[WGS_1984,
SPHEROID[WGS 84,6378137,298.2572235629972,
AUTHORITY[EPSG,7030]],
AUTHORITY[EPSG,6326]],
PRIMEM[Greenwich,0],
UNIT[degree,0.0174532925199433],
AUTHORITY[EPSG,4326]]
Origin = (65.000,5.0028436)
Pixel Size = (0.0008333,0.0008333)
Metadata: snip

Thanks for confirming that there is a bug in reading coordinate system
info from netCDF files.  We're trying to make sure that some of our
netCDF's in Lambert projection are correctly written - specifically so
I can read them with gdal, and could not figure out what values to
assign them.

I think you're very close to having the netCDF driver working well
enough to be used for serving and consuming WCS, but there are still a
few issues to work out.  What can I do to help?  Setup sample
problems, test, write use-cases?  I'm trying  to get a MapServer-based
WCS system up at work that uses tiled netCDF's with time series data,
so I think I've hit most of the remaining issues along the way.

Roger
--

On Mon, Mar 16, 2009 at 9:51 AM, Frank Warmerdam warmer...@pobox.com wrote:
 Roger André wrote:

 Hi all,

 Just a couple quick netCDF related questions:

 1)  Is there any way to include projection information in a netCDF
 file so that gdal can properly read it?  It appears that
 georeferencing is handled correctly, but not projection info.  This
 causes problems with the whole north/south origin shift between the
 netCDF and other raster formats when gdal is used to convert the .nc
 file into a GeoTIFF.  Most notably, I'm not sure how to assign the
 data a projection - either for display in MapServer, or to reproject
 it to EPSG:4326:  See example below.

 File created with ArcGIS 9.3:
 $ gdalinfo fwdir_bas.nc
 Driver: netCDF/Network Common Data Format
 Files: fwdir_bas.nc
 Size is 42000, 42000
 Coordinate System is `'
 Origin = (65.000,5.0028436)
 Pixel Size = (0.0008333,0.0008333)
 Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#Source_Software=ESRI ArcGIS
  Fwdir_bas#long_name=Fwdir_bas

  Fwdir_bas#esri_pe_string=GEOGCS[GCS_WGS_1984,DATUM[D_WGS_1984,SPHEROID[WGS_1984,6378137.0,298.257223563]],PRIMEM[Greenwich,0.0],UNIT[Degree,0.0174532925199433]]
  Fwdir_bas#coordinates=lon lat
  Fwdir_bas#units=Degree
  Fwdir_bas#missing_value=
  lon#long_name=longitude coordinate
  lon#standard_name=longitude
  lon#units=degrees_east
  lat#long_name=latitude coordinate
  lat#standard_name=latitude
  lat#units=degrees_north

 Roger,

 The above example is not projected as far as I can see.  It is
 geographic.  So I'm not sure what this was intended to be an
 example of.

 If I use gdal_translate to convert a utm GeoTIFF file to
 Netcdf I get a gdalinfo report like below.

 Driver: netCDF/Network Common Data Format
 Files: out.nc
 Size is 512, 512
 Coordinate System is:
 PROJCS[NAD27 / UTM zone 11N,
    GEOGCS[NAD27,
        DATUM[North_American_Datum_1927,
            SPHEROID[Clarke 1866,6378206.4,294.9786982139006,
                AUTHORITY[EPSG,7008]],
            AUTHORITY[EPSG,6267]],
        PRIMEM[Greenwich,0],
        UNIT[degree,0.0174532925199433],
        AUTHORITY[EPSG,4267]],
    PROJECTION[Transverse_Mercator],
    PARAMETER[latitude_of_origin,0],
    PARAMETER[central_meridian,-117],
    PARAMETER[scale_factor,0.9996],
    PARAMETER[false_easting,50],
    PARAMETER[false_northing,0],
    UNIT[metre,1,
        AUTHORITY[EPSG,9001]],
    AUTHORITY[EPSG,26711]]
 Origin = (440720.000,3751320.000)
 Pixel Size = (60.000,-60.000)
 Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#AREA_OR_POINT=Area
  Band1#grid_mapping=transverse_mercator
  Band1#long_name=GDAL Band Number 1
  transverse_mercator#Northernmost_Northing=3.75132e+06
  transverse_mercator#Southernmost_Northing=3.7206e+06
  transverse_mercator#Easternmost_Easting=471440
  transverse_mercator#Westernmost_Easting=440720
  transverse_mercator#spatial_ref=PROJCS[NAD27 / UTM zone
 11N,GEOGCS[NAD27,DATUM[North_American_Datum_1927,SPHEROID[Clarke
 1866,6378206.4,294.9786982139006,AUTHORITY[EPSG,7008]],AUTHORITY[EPSG,6267]],PRIMEM[Greenwich,0],UNIT[degree,0.0174532925199433],AUTHORITY[EPSG,4267]],PROJECTION[Transverse_Mercator],PARAMETER[latitude_of_origin,0],PARAMETER[central_meridian,-117],PARAMETER[scale_factor,0.9996],PARAMETER[false_easting

[gdal-dev] Re: gdal2tiles dateline gap

[Roger André]

Hi Jeff,

Adding list, just in case there is additional input on this.

Yes, but the solution I found may not be applicable to your use-case.
I needed to create Google Map tiles in spherical mercator projection,
not Google Earth KML, so I was able to use TileCache with MapServer to
generate them.  I found that I had more flexibility using that
solution, than I did using the gdal2tiles tool.

I should also add that when I used input data that was already in
sperical mercator projection, I was unable to get data all the way out
to the edges of their defined extents.  You might try converting the
input data to EPSG:4326 first, then see what you get out of
gdal2tiles.  I believe it can handle the reprojection for you
on-the-fly.

Roger
--

On Sun, Mar 15, 2009 at 4:41 PM, Jeff Logan jlo...@customweather.com wrote:
 Hi Roger,

 Gdal2tiles output  is chopping off New Zealand from a global data set...the
 input spherical mercator global image has New Zealand intact. Did you
 discover a solution after posting to 'gdal.dev' nine weeks ago?

 Thanks, Jeff


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

[gdal-dev] gdaltindex and netCDF files

[Roger André]

Hi All,

Is it possible to use gdaltindex to create a shapefile index of a
group of netCDF files?  Specifivally, I need to index a specific
netCDF variable that exists in all of the files.  I just tried running
it on a single file, using the variable selection syntax that seems to
work with the other gdal utilities, but can't seem to make it work.

gdaltindex netcdf_index.shp
./NETCDF:N35E115_5x5_1997010100_2007010223.ttnc:elev

ERROR 4: `./NETCDF:N35E115_5x5_1997010100_2007010223.ttnc:elev'
does not exist in the file system,
and is not recognised as a supported dataset name.

Unable to open ./NETCDF:N35E115_5x5_1997010100_2007010223.ttnc:elev,
skipping.

Thanks,

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

[gdal-dev] MapServer using GDAL_WMS to access TileCache layer

[Roger André]

Hi All,

I've created a TileCache layer that's being used to overlay tiles in a
Google Maps page.  That's working just fine, but I would also like to
access this layer via a MapServer WMS call as well.  I have tried to
setup a GDAL_WMS file that can read the TileCache layer, but am
getting IReadBlock failed errors when I try to run gdal_translate on
it, or try to view the layer in a MapServer map.

Below are the pertinent details of my configuration.

TileCache layer:
---
[20m_global]
type=MapServerLayer
mapfile=/var/www/mapfiles/20m_global/20m_global.map
#url=http://10.10.40.152/cgi-bin/mapserv?map=/var/www/mapfiles/20m_global/20m_global.map
metaTile=true
metaSize=2,2
metaBuffer=0,500
layers=imagery
spherical_mercator=true
srs=EPSG:900913

GDAL_WMS xml file:
--
GDAL_WMS
Service name=TileService
!--ServerUrlhttp://tile.openstreetmap.org/${z}/${x}/${y}.png/ServerUrl--
ServerUrlhttp://10.10.10.143/cgi-bin/tilecache-2.04/tilecache.cgi/1.0.0/v2_80m_global/{Z}/{X}/{Y}.png/ServerUrl
/Service
DataWindow
UpperLeftX-20037508.34/UpperLeftX
UpperLeftY20037508.34/UpperLeftY
LowerRightX20037508.34/LowerRightX
LowerRightY-20037508.34/LowerRightY
TileLevel11/TileLevel
TileCountX1/TileCountX
TileCountY1/TileCountY
YOrigintop/YOrigin
/DataWindow
ProjectionEPSG:900913/Projection
BlockSizeX256/BlockSizeX
BlockSizeY256/BlockSizeY
BandsCount3/BandsCount
Cache/
/GDAL_WMS

gdalinfo output:
-
$ gdalinfo tilecache.xml
Driver: WMS/OGC Web Map Service
Files: tilecache.xml
Size is 524288, 524288
Coordinate System is:
PROJCS[Google Maps Global Mercator,
GEOGCS[WGS 84,
DATUM[WGS_1984,
SPHEROID[WGS 84,6378137,298.257223563,
AUTHORITY[EPSG,7030]],
AUTHORITY[EPSG,6326]],
PRIMEM[Greenwich,0,
AUTHORITY[EPSG,8901]],
UNIT[degree,0.01745329251994328,
AUTHORITY[EPSG,9122]],
AUTHORITY[EPSG,4326]],
PROJECTION[Mercator_2SP],
PARAMETER[standard_parallel_1,0],
PARAMETER[latitude_of_origin,0],
PARAMETER[central_meridian,0],
PARAMETER[false_easting,0],
PARAMETER[false_northing,0],
UNIT[Meter,1],
EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgri...@null +wktext
+no_defs],
AUTHORITY[EPSG,900913]]
Origin = (-20037508.33999850988,20037508.33999850988)
Pixel Size = (76.437028274536132,-76.437028274536132)
Corner Coordinates:
Upper Left  (-20037508.340,20037508.340) (180d 0'0.00W, 85d 3'4.06N)
Lower Left  (-20037508.340,-20037508.340) (180d 0'0.00W, 85d 3'4.06S)
Upper Right (20037508.340,20037508.340) (180d 0'0.00E, 85d 3'4.06N)
Lower Right (20037508.340,-20037508.340) (180d 0'0.00E, 85d 3'4.06S)
Center  (   0.000,   0.000) (  0d 0'0.01E,  0d 0'0.00S)
Band 1 Block=256x256 Type=Byte, ColorInterp=Undefined
  Overviews: 262144x262144, 131072x131072, 65536x65536, 32768x32768,
16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512,
256x256
Band 2 Block=256x256 Type=Byte, ColorInterp=Undefined
  Overviews: 262144x262144, 131072x131072, 65536x65536, 32768x32768,
16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512,
256x256
Band 3 Block=256x256 Type=Byte, ColorInterp=Undefined
  Overviews: 262144x262144, 131072x131072, 65536x65536, 32768x32768,
16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512,
256x256

gdal_translate error:

$ gdal_translate -outsize 1% 1% tilecache.xml test.tif
Input file size is 524288, 524288
0ERROR 4: `/vsimem/wms/0x8503470/wmsresult.dat' not recognised as a
supported file format.

ERROR 1: GDALWMS: Unable to open downloaded block.
ERROR 1: GDALWMS: ReadBlockFromFile failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 0
ERROR 4: `/vsimem/wms/0x809dc08/wmsresult.dat' not recognised as a
supported file format.

ERROR 1: GDALWMS: Unable to open downloaded block.
ERROR 1: GDALWMS: ReadBlockFromFile failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 0

MapServer error:
---
msDrawMap(): Image handling error. Failed to draw layer named
'imagery'. IReadBlock failed at X offset 1442, Y offset 607: Unable to
access file. GDALRasterIO() failed: drawGDAL()

If you see anything that jumps out at you as being incorrect, I'd love
to hear from you.

Thanks,

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

Re: [gdal-dev] MapServer using GDAL_WMS to access TileCache layer

[Roger André]

Hi Even,

Thanks for the info.  I have updated my GDAL_WMS file to now show:

ServerUrlhttp://10.10.10.143/cgi-bin/tilecache-2.04/tilecache.cgi/1.0.0/v2_80m_global/${z}/${x}/${y}.png/ServerUrl


Still doesn't work though.


gdal_translate -outsize 1% 1% tilecache.xml test.tif
Input file size is 524288, 524288
0ERROR 4: `/vsimem/wms/0x8503480/wmsresult.dat' not recognised as a
supported file format.

ERROR 1: GDALWMS: Unable to open downloaded block.
ERROR 1: GDALWMS: ReadBlockFromFile failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 0
ERROR 4: `/vsimem/wms/0x809dc28/wmsresult.dat' not recognised as a
supported file format.

ERROR 1: GDALWMS: Unable to open downloaded block.
ERROR 1: GDALWMS: ReadBlockFromFile failed.
ERROR 1: IReadBlock failed at X offset 0, Y offset 0


On Mon, Mar 23, 2009 at 1:35 PM, Even Rouault
even.roua...@mines-paris.org wrote:
 Roger,

 your GDAL_WMS xml file is not correct. The format of the URL you use is not
 consistant with the one expected by the TileService service. You've used
 the example of TMS service (which really new in GDAL 1.7.0dev SVN only) but
 the correct format would be exactly what you have commented out (${x} in
 lowercase, and not {X}).

 Best regards,
 Even

 Le Monday 23 March 2009 21:18:16 Roger André, vous avez écrit :
 Hi All,

 I've created a TileCache layer that's being used to overlay tiles in a
 Google Maps page.  That's working just fine, but I would also like to
 access this layer via a MapServer WMS call as well.  I have tried to
 setup a GDAL_WMS file that can read the TileCache layer, but am
 getting IReadBlock failed errors when I try to run gdal_translate on
 it, or try to view the layer in a MapServer map.

 Below are the pertinent details of my configuration.

 TileCache layer:
 ---
 [20m_global]
 type=MapServerLayer
 mapfile=/var/www/mapfiles/20m_global/20m_global.map
 #url=http://10.10.40.152/cgi-bin/mapserv?map=/var/www/mapfiles/20m_global/2
0m_global.map metaTile=true
 metaSize=2,2
 metaBuffer=0,500
 layers=imagery
 spherical_mercator=true
 srs=EPSG:900913

 GDAL_WMS xml file:
 --
 GDAL_WMS
 Service name=TileService
 !--ServerUrlhttp://tile.openstreetmap.org/${z}/${x}/${y}.png/ServerUrl
--
 ServerUrlhttp://10.10.10.143/cgi-bin/tilecache-2.04/tilecache.cgi/1.0.0/v
2_80m_global/{Z}/{X}/{Y}.png/ServerUrl /Service
 DataWindow
 UpperLeftX-20037508.34/UpperLeftX
 UpperLeftY20037508.34/UpperLeftY
 LowerRightX20037508.34/LowerRightX
 LowerRightY-20037508.34/LowerRightY
 TileLevel11/TileLevel
 TileCountX1/TileCountX
 TileCountY1/TileCountY
 YOrigintop/YOrigin
 /DataWindow
 ProjectionEPSG:900913/Projection
 BlockSizeX256/BlockSizeX
 BlockSizeY256/BlockSizeY
 BandsCount3/BandsCount
 Cache/
 /GDAL_WMS

 gdalinfo output:
 -
 $ gdalinfo tilecache.xml
 Driver: WMS/OGC Web Map Service
 Files: tilecache.xml
 Size is 524288, 524288
 Coordinate System is:
 PROJCS[Google Maps Global Mercator,
     GEOGCS[WGS 84,
         DATUM[WGS_1984,
             SPHEROID[WGS 84,6378137,298.257223563,
                 AUTHORITY[EPSG,7030]],
             AUTHORITY[EPSG,6326]],
         PRIMEM[Greenwich,0,
             AUTHORITY[EPSG,8901]],
         UNIT[degree,0.01745329251994328,
             AUTHORITY[EPSG,9122]],
         AUTHORITY[EPSG,4326]],
     PROJECTION[Mercator_2SP],
     PARAMETER[standard_parallel_1,0],
     PARAMETER[latitude_of_origin,0],
     PARAMETER[central_meridian,0],
     PARAMETER[false_easting,0],
     PARAMETER[false_northing,0],
     UNIT[Meter,1],
     EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgri...@null +wktext
 +no_defs],
     AUTHORITY[EPSG,900913]]
 Origin = (-20037508.33999850988,20037508.33999850988)
 Pixel Size = (76.437028274536132,-76.437028274536132)
 Corner Coordinates:
 Upper Left  (-20037508.340,20037508.340) (180d 0'0.00W, 85d 3'4.06N)
 Lower Left  (-20037508.340,-20037508.340) (180d 0'0.00W, 85d 3'4.06S)
 Upper Right (20037508.340,20037508.340) (180d 0'0.00E, 85d 3'4.06N)
 Lower Right (20037508.340,-20037508.340) (180d 0'0.00E, 85d 3'4.06S)
 Center      (   0.000,   0.000) (  0d 0'0.01E,  0d 0'0.00S)
 Band 1 Block=256x256 Type=Byte, ColorInterp=Undefined
   Overviews: 262144x262144, 131072x131072, 65536x65536, 32768x32768,
 16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512,
 256x256
 Band 2 Block=256x256 Type=Byte, ColorInterp=Undefined
   Overviews: 262144x262144, 131072x131072, 65536x65536, 32768x32768,
 16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512,
 256x256
 Band 3 Block=256x256 Type=Byte, ColorInterp=Undefined
   Overviews: 262144x262144, 131072x131072, 65536x65536, 32768x32768,
 16384x16384, 8192x8192, 4096x4096, 2048x2048, 1024x1024, 512x512,
 256x256

 gdal_translate error:
 
 $ gdal_translate -outsize 1% 1% tilecache.xml test.tif
 Input file size is 524288, 524288
 0ERROR 4: `/vsimem/wms/0x8503470

[gdal-dev] GDAL TMS minidriver

[Roger André]

Hi All,

I'm trying to use the TMS minidriver to read a TileCache instance and am
getting an error that says, GDALWMS: No mini-driver registered for 'TMS'.
I'm trying to follow the instructions on http://www.gdal.org/frmt_wms.html,
but am a bit confused. At the top of the page it says there are 3 services
which are currently supported; WMS, WorldWind, and TileService.  However,
lower in the doc it also mentions TMS, and one of the examples
http://www.gdal.org/frmt_wms_openstreetmap_tms.xml uses that service type,
and looks like it would work with TileCache.

I'm using gdal 1.5.2, is this functionality supported in that version?

Thanks,

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

Re: [gdal-dev] GDAL TMS minidriver

[Roger André]

Got it, V 1.7.

Thanks.
--

2009/3/26 Klokan Petr Přidal klo...@klokan.cz

 Roger,

 this patch was done by Christopher Schmidt during Toronto Code Sprint
 2009 in the beginning of this month.
 The functionality is available only in gdal-trunk (1.7dev) since r16585.

 Details:
 http://trac.osgeo.org/gdal/ticket/2878
 http://crschmidt.net/blog/376/toronto-code-sprint-day-2/

 If you find a bug in the online documentation feel free to patch it ;-).

 Best reagards,

 Klokan

 2009/3/26 Roger André ran...@gmail.com:
  Hi All,
 
  I'm trying to use the TMS minidriver to read a TileCache instance and am
  getting an error that says, GDALWMS: No mini-driver registered for
 'TMS'.
  I'm trying to follow the instructions on
 http://www.gdal.org/frmt_wms.html,
  but am a bit confused. At the top of the page it says there are 3
 services
  which are currently supported; WMS, WorldWind, and TileService.  However,
  lower in the doc it also mentions TMS, and one of the examples
  http://www.gdal.org/frmt_wms_openstreetmap_tms.xml uses that service
 type,
  and looks like it would work with TileCache.
 
  I'm using gdal 1.5.2, is this functionality supported in that version?
 
  Thanks,
 
  Roger
  --
 
  ___
  gdal-dev mailing list
  gdal-dev@lists.osgeo.org
  http://lists.osgeo.org/mailman/listinfo/gdal-dev
 



 --
 http://blog.klokan.cz/
 http://www.maptiler.org/
 http://www.oldmapsonline.org/

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

Re: [gdal-dev] Problem in importing to Google Earth a nitf image converted from a tif image by gdal_translate

[Roger André]

Hi Mo,

Thanks for that pointer on the differences between GE Layers and Places.
I'll be sure to take a close look at those links.  I've had occasion in the
past to view georeferenced imagery in Google Earth, much of it created or
altered by gdal in some way, so this topic interests me.

I'm curious, is your goal to view the NITF format specifically and to have
GE read the internal georeferencing info, or is it just to view the image in
GE?  It sounds like you have already accomplished the latter using a KML
Ground Overlay and a tif, a method which has also worked reliably for me in
the past as well. I've also had success with PNG and JPEG format files.

Since in your .nitf test the georeferencing was supplied by the KML file,
the red X you are getting implies to me that the image format isn't being
correctly interpreted by GE for some reason.  The fact that Envi chokes on
it as well does seem to point to gdal as a possible culprit, although the
NITF format is not one that I have seen very much.  Have you tried importing
the GeoTIFF version of your file directly into GE Pro as well?  Also, do you
have a .nitf file that was not created by gdal, and which you could use to
test GE via the Add Image Overlay tool?

Roger
--

On Thu, Apr 9, 2009 at 8:37 AM, mohwawang mohwaw...@yahoo.com wrote:


 Roger,

 I'm using Google Earth Pro 5.0. Thank you for the suggestions.

 Google Earth Pro 4.0 also supports importing NITF:
 http://earth.google.com/userguide/v4/ug_importdata.html that is for
 importing your data (including vector data and imagery data) to the Places
 of Google Earth while
 http://earth.google.com/userguide/v4/ug_imageoverlays.html is for Using
 Image Overlays and 3D Models.

 See http://earth.google.com/userguide/v4/ug_layers.html#layers_vs_placesfor 
 the difference between Layers and Places in GE.

 This is the first time I'm importing geo-referenced images to Google Earth
 Pro and so far I don't have much luck. As I said that I could not use GE to
 import the NITF images converted from gdal_translate (with EPSG:4326 that is
 the same projection, Simple Cylindrical projection, or known as Lat/Lon
 WGS84 that GE uses). So I tried to narrow down the problem. This time I
 imported to GE an native Erdas IMAGINE image (.img) (with UTM projection
 that GE supports too) but still failed. The later test case has nothing to
 do with gdal and is purely a question for GE. That .img image can be opened
 by both IMAGINE and OpenEV.

 The NITF image (converted by gdal_translate) not just could not be opened
 by GE it could not be opened by ENVI 4.1 either. It can be opened by OpenEV
 (tightly related to gdal) though. So it could be a problem in either
 gdal_translate or ENVI (I doubt though) or both.

 Thanks.

 Mo

 --- On Thu, 4/9/09, Roger André ran...@gmail.com wrote:

  From: Roger André ran...@gmail.com
  Subject: Re: [gdal-dev] Problem in importing to Google Earth a nitf image
  converted from a tif image by gdal_translate
  To: mohwawang mohwaw...@yahoo.com
  Cc: gdal-dev@lists.osgeo.org gdal-dev@lists.osgeo.org
  Date: Thursday, April 9, 2009, 12:44 AM
  Granted, I haven't tried V 5 yet,
  but your quote also states, Google Earth Pro or
  EC.  Are you using one of those versions of the
  client?  In any case, if the file opens in OpenEV, chances
  are good that the .ntf generated by gdal is ok, and that
  there is some sort of problem with how Google Earth is
  reading it.  Have you tried emailing the KML support list
  as well?
 
  --
 
  On Wed, Apr 8, 2009 at 8:51 PM,
  mohwawang mohwaw...@yahoo.com
  wrote:
 
 
 
  Google Earth 5.0 supports NITF see
 
  http://earth.google.com/intl/en/userguide/v5/ug_importdata.html
 
 
 
  [Importing Imagery
 
 
 
  You can open GIS imagery files in Google Earth Pro or EC
  and have the files correctly projected over the proper map
  coordinates in the 3D viewer. Google Earth Pro and EC
  support the following types of GIS imagery:
 
 
 
  * TIFF (.tif), including GeoTiff and compressed TIFF
  files
 
  * National Imagery Transmission Format (.ntf)
 
  * Erdas Imagine Images (.img)
 
 
 
  In addition, you can also import the following images.
  Because they do not contain projection information, you must
  manually edit their coordinates for correct positioning:
 
 
 
  * Atlantis MFF Raster (.hdr)
 
  * PCIDSK Database File (.pix)
 
  * Portable Pixmap Format (.pnm)
 
  * Device Independent Bitmap (.bmp)
 
 
 
  Note - All imagery files must contain the correct
  projection information in order to be accurately
  re-projected by the Google Earth software.]
 
 
 
 
 
  Mo
 
 
 
  --- On Wed, 4/8/09, Roger André ran...@gmail.com
  wrote:
 
 
 
   From: Roger André ran...@gmail.com
 
   Subject: Re: [gdal-dev] Problem in importing to Google
  Earth a nitf image  converted from a tif image by
  gdal_translate
 
   To: mohwawang mohwaw...@yahoo.com
 
   Cc: gdal-dev@lists.osgeo.org
  gdal-dev@lists.osgeo.org
 
   Date: Wednesday, April 8

Re: [gdal-dev] Problem in importing to Google Earth a nitf image converted from a tif image by gdal_translate

[Roger André]

Hi Mo,

What is the gdalinfo output for the GeoTIFF that you created?  If you didn't
get an error when you tried to load it in GE, perhaps the georeferencing in
the image doesn't quite match what you specified in your KML.

Roger
--

On Thu, Apr 9, 2009 at 1:21 PM, mohwawang mohwaw...@yahoo.com wrote:


 Hi Roger,

 My goal is to be able to directly import NITF image to GE without using
 KML. You can see the advantage of this application - if it works. I'm glad
 that you have viewed in GE georeferenced imagery created or altered by gdal.
 Maybe you can point out what I did wrong.

 It is good that you asked about GeoTIFF. So I used gdal_translate to
 convert the tif image to a GeoTIFF image. I tried to drag and drop the
 GeoTIFF image to GE but nothing happened. I did not see the GeoTIFF image in
 GE and I did not get any error message either.

 I can open the GeoTIFF image with OpenEV and I can read the lan/lon
 coordinates and they are correct. One thing is that it seems that OpenEV
 does not set North to Up. Its orientation aligns with the file coordinates -
 (0,0) is at the upper left corner. I then opened tif file in the same OpenEV
 and they match to each other.

 Again, I can import the tif image with its kml to GE but I can not import
 the GeoTIFF to GE (without using kml), it seems that nothing happened. I can
 import the GeoTIFF to GE with kml - though that is not what I want. But the
 GeoTIFF is better than the NITF that could not be imported to GE even with
 the same kml.

 What's wrong with GeoTIFF? I mean why I could not import it to GE? I listed
 the gdal_translate command line here:

 gdal_translate -of GTiff tiff.tif gtiff.gtf -a_srs EPSG:4326 -gcp 0 0
 -70.887525562763 42.8093748999134 -gcp 560 0 -70.887525562763
 42.8082925732318 -gcp 560 419 -70.8894939447889 42.8082925732318 -gcp 0 419
 -70.8894939447889 42.8093748999134

 I don't have a .ntif file that was not created by gdal to test GE via the
 Import Image Tool.

 Thanks.

 Mo


 --- On Thu, 4/9/09, Roger André ran...@gmail.com wrote:

  From: Roger André ran...@gmail.com
  Subject: Re: [gdal-dev] Problem in importing to Google Earth a nitf image
  converted from a tif image by gdal_translate
  To: mohwawang mohwaw...@yahoo.com
  Cc: gdal-dev@lists.osgeo.org gdal-dev@lists.osgeo.org
  Date: Thursday, April 9, 2009, 12:47 PM
  Hi Mo,
 
  Thanks for that pointer on the differences between GE
  Layers and Places.  I'll be sure to take a close look
  at those links.  I've had occasion in the past to view
  georeferenced imagery in Google Earth, much of it created or
  altered by gdal in some way, so this topic interests me.
 
 
  I'm curious, is your goal to view the NITF format
  specifically and to have GE read the internal georeferencing
  info, or is it just to view the image in GE?  It sounds
  like you have already accomplished the latter using a KML
  Ground Overlay and a tif, a method which has also worked
  reliably for me in the past as well. I've also had
  success with PNG and JPEG format files.
 
 
  Since in your .nitf test the georeferencing was supplied by
  the KML file, the red X you are getting implies
  to me that the image format isn't being correctly
  interpreted by GE for some reason.  The fact that Envi
  chokes on it as well does seem to point to gdal as a
  possible culprit, although the NITF format is not one that I
  have seen very much.  Have you tried importing the GeoTIFF
  version of your file directly into GE Pro as well?  Also,
  do you have a .nitf file that was not created by gdal, and
  which you could use to test GE via the Add Image
  Overlay tool?
 
 
  Roger
  --
 
  On Thu, Apr 9, 2009 at 8:37 AM,
  mohwawang mohwaw...@yahoo.com
  wrote:
 
 
 
  Roger,
 
 
 
  I'm using Google Earth Pro 5.0. Thank you for the
  suggestions.
 
 
 
  Google Earth Pro 4.0 also supports importing NITF:
 
  http://earth.google.com/userguide/v4/ug_importdata.html
  that is for importing your data (including vector data and
  imagery data) to the Places of Google Earth while
 http://earth.google.com/userguide/v4/ug_imageoverlays.html
  is for Using Image Overlays and 3D Models.
 
 
 
 
  See http://earth.google.com/userguide/v4/ug_layers.html#layers_vs_places
  for the difference between Layers and Places in GE.
 
 
 
 
  This is the first time I'm importing geo-referenced
  images to Google Earth Pro and so far I don't have much
  luck. As I said that I could not use GE to import the NITF
  images converted from gdal_translate (with EPSG:4326 that is
  the same projection, Simple Cylindrical projection, or known
  as Lat/Lon WGS84 that GE uses). So I tried to narrow down
  the problem. This time I imported to GE an native Erdas
  IMAGINE image (.img) (with UTM projection that GE supports
  too) but still failed. The later test case has nothing to do
  with gdal and is purely a question for GE. That .img image
  can be opened by both IMAGINE and OpenEV.
 
 
 
 
  The NITF image (converted by gdal_translate) not just

Re: [Gdal-dev] 1-band TIFF FILE color ramp NEEDED

[Roger André]

Here's another way you can do it, a Python script that will create discreet
color classes from a grayscale image.  It uses a function named MakeColor
which is currently set to work on Z values between 0 and 255.  You can test
it by converting your file into an 8-bit version (gdal_translate -scale -ot
Byte infile outfile ), and then running the tool on it.  Once you see
how it works, it should be pretty simple to modify the function to deal with
the actual ranges in your data.

Roger
--

On Thu, Apr 9, 2009 at 4:24 PM, alvarez00 oalvare...@gmail.com wrote:


 I'm trying to make a color ramp for a 1-band tiff file or a jpeg. I have a
 1-band tiff file and I converted that tiff file into VRT. I was reading some
 of the posts and some suggest to modify the VRT file and add a ColorTable
 tag but that doesn't work. My objective is to show the tiff(or)jpeg on
 google earth. I'm running version 1.4.1.0 of gdal. Below is the VRT file
 without any modifications. What I modify was the
 'ColorInterpGray/ColorInterp' to 'ColorInterpPalette/ColorInterp'
 and added a ColorTable.

 VRTDataset rasterXSize=1200 rasterYSize=1096
  
 SRSGEOGCS[quot;NAD83quot;,DATUM[quot;North_American_Datum_1983quot;,SPHEROID[quot;GRS
 1980quot;,6378137,298.2572221010002,AUTHORITY[quot;EPSGquot;,quot;7019quot;]],AUTHORITY[quot;EPSGquot;,quot;6269quot;]],PRIMEM[quot;Greenwichquot;,0],UNIT[quot;degreequot;,0.0174532925199433],AUTHORITY[quot;EPSGquot;,quot;4269quot;]]/SRS
  GeoTransform -1.2459079918952935e+02,  8.9964705882352915e-03,
  0.e+00,  4.2216820949146452e+01,  0.e+00,
 -8.9964705882352915e-03/GeoTransform
  Metadata
MDI key=AREA_OR_POINTArea/MDI
  /Metadata
  VRTRasterBand dataType=Float32 band=1
Metadata/
NoDataValue-3.40282346638529E+38/NoDataValue
ColorInterpGray/ColorInterp
SimpleSource
  SourceFilename
 relativeToVRT=1GoesWest1V1561915.tif/SourceFilename
  SourceBand1/SourceBand
  SrcRect xOff=0 yOff=0 xSize=1200 ySize=1096/
  DstRect xOff=0 yOff=0 xSize=1200 ySize=1096/
/SimpleSource
  /VRTRasterBand
 /VRTDataset

 --
 View this message in context:
 http://n2.nabble.com/1-band-TIFF-FILE-color-ramp-NEEDED-tp2613963p2613963.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

#! /usr/bin/env python

# Create a colored image based on LUT in MakeColor function

# Author: Roger Andre, October 2008

# Usage: discreet_gray2color.py infile outfile

from osgeo import gdal
import sys
import numpy
import os.path
gdal.TermProgress = gdal.TermProgress_nocb

src_file = sys.argv[1]
dst_file = sys.argv[2]
out_bands = 3

def MakeColor(z_value):
  '''LUT for color ramp. Keys are pixel values, hash values are RGB triplets.'''
  color_dict = {
0:[102,0,255],
   25:[20,82,255],
   50:[0,194,224],
   75:[0,255,122],
  100:[41,255,0],
  125:[204,255,0],
  150:[245,194,0],
  175:[224,118,0],
  200:[168,46,0],
  225:[105,0,0],
  250:[64,0,0],
  255:[0,0,0]}

  key_list = color_dict.keys()
  key_list.sort()
  while len(key_list)  0:
last_val = key_list[-1]
if z_value = last_val:
  return color_dict[last_val]
else:
  key_list.remove(last_val)

# Print some info
print Creating %s % (dst_file)

# Open source file
src_ds = gdal.Open( src_file )
src_band = src_ds.GetRasterBand(1)

# create destination file
dst_driver = gdal.GetDriverByName('GTiff')
dst_ds = dst_driver.Create(dst_file, src_ds.RasterXSize, src_ds.RasterYSize, out_bands, gdal.GDT_Byte) 

# create output bands
band1 = numpy.zeros([src_ds.RasterYSize, src_ds.RasterXSize])
band2 = numpy.zeros([src_ds.RasterYSize, src_ds.RasterXSize])
band3 = numpy.zeros([src_ds.RasterYSize, src_ds.RasterXSize])

# set the projection and georeferencing info
dst_ds.SetProjection( src_ds.GetProjection() )
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )

# read the source file
gdal.TermProgress( 0.0 )
for iY in range(src_ds.RasterYSize):
  src_data = src_band.ReadAsArray(0,iY,src_ds.RasterXSize,1)
  col_values = src_data[0] # array of z_values, one per row in source data
  for iX in range(src_ds.RasterXSize):
z_value = col_values[iX]
# print z_value # uncomment to see what value breaks color ramp
[R,G,B] = MakeColor(z_value)
band1[iY][iX] = R
band2[iY][iX] = G
band3[iY][iX] = B
  gdal.TermProgress( (iY+1.0) / src_ds.RasterYSize )

# write each band out
dst_ds.GetRasterBand(1).WriteArray(band1)
dst_ds.GetRasterBand(2).WriteArray(band2)
dst_ds.GetRasterBand(3).WriteArray(band3)

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

Re: [Gdal-dev] 1-band TIFF FILE color ramp NEEDED

[Roger André]

Sorry about that.  It has to do with the little progress indicator that runs
when gdal processes a file.  There are different ways of invoking it, and my
mileage varies with each of them.  I think it is based on what version of
gdal you're using.  Anyhow, 2 ways you can deal with this.  I think with
your version of gdal you can just remove the import statement, line 14 (
gdal.TermProgress = gdal.TermProgress_nocb ), entirely. If that doesn't
work, just remove all 3 lines that have TermProgress in them.  You won't
have a status indicator then to tell you how far along the coloring is, but
the program should run.  Also, just fyi that the colorization won't be
incredibly fast. It takes around 5 minutes to color an image that is 2500 x
2500 pixels on my machine.  There is obviously room for optimization here
;).
--

On Sat, Apr 11, 2009 at 5:05 PM, alvarez00 oalvare...@gmail.com wrote:


 Hi Roger,

 I'm trying to run the python code you suggested but I keep getting this
 error :
 Traceback (most recent call last):
  File discreet_gray2color.py, line 14, in module
gdal.TermProgress = gdal.TermProgress_nocb
 AttributeError: 'module' object has no attribute 'TermProgress_nocb' 

 I don't know if it's the version of gdal I'm using or something else. (i'm
 using gdal 1.4.1.0)

 thanks,




 Here's another way you can do it, a Python script that will create discreet
 color classes from a grayscale image.  It uses a function named MakeColor
 which is currently set to work on Z values between 0 and 255.  You can test
 it by converting your file into an 8-bit version (gdal_translate -scale -ot
 Byte infile outfile ), and then running the tool on it.  Once you see
 how it works, it should be pretty simple to modify the function to deal
 with
 the actual ranges in your data.

 Roger
 --

 On Thu, Apr 9, 2009 at 4:24 PM, alvarez00 oalvare...@gmail.com wrote:

 
  I'm trying to make a color ramp for a 1-band tiff file or a jpeg. I have
 a
  1-band tiff file and I converted that tiff file into VRT. I was reading
 some
  of the posts and some suggest to modify the VRT file and add a ColorTable
  tag but that doesn't work. My objective is to show the tiff(or)jpeg on
  google earth. I'm running version 1.4.1.0 of gdal. Below is the VRT file
  without any modifications. What I modify was the
  'ColorInterpGray/ColorInterp' to 'ColorInterpPalette/ColorInterp'
  and added a ColorTable.
 
  VRTDataset rasterXSize=1200 rasterYSize=1096
 
  
 SRSGEOGCS[quot;NAD83quot;,DATUM[quot;North_American_Datum_1983quot;,SPHEROID[quot;GRS
 
 1980quot;,6378137,298.2572221010002,AUTHORITY[quot;EPSGquot;,quot;7019quot;]],AUTHORITY[quot;EPSGquot;,quot;6269quot;]],PRIMEM[quot;Greenwichquot;,0],UNIT[quot;degreequot;,0.0174532925199433],AUTHORITY[quot;EPSGquot;,quot;4269quot;]]/SRS
   GeoTransform -1.2459079918952935e+02,  8.9964705882352915e-03,
   0.e+00,  4.2216820949146452e+01,
  0.e+00,
  -8.9964705882352915e-03/GeoTransform
   Metadata
 MDI key=AREA_OR_POINTArea/MDI
   /Metadata
   VRTRasterBand dataType=Float32 band=1
 Metadata/
 NoDataValue-3.40282346638529E+38/NoDataValue
 ColorInterpGray/ColorInterp
 SimpleSource
   SourceFilename
  relativeToVRT=1GoesWest1V1561915.tif/SourceFilename
   SourceBand1/SourceBand
   SrcRect xOff=0 yOff=0 xSize=1200 ySize=1096/
   DstRect xOff=0 yOff=0 xSize=1200 ySize=1096/
 /SimpleSource
   /VRTRasterBand
  /VRTDataset
 
  --
  View this message in context:
 
 http://n2.nabble.com/1-band-TIFF-FILE-color-ramp-NEEDED-tp2613963p2613963.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
 

 #! /usr/bin/env python

 # Create a colored image based on LUT in MakeColor function

 # Author: Roger Andre, October 2008

 # Usage: discreet_gray2color.py infile outfile

 from osgeo import gdal
 import sys
 import numpy
 import os.path
 gdal.TermProgress = gdal.TermProgress_nocb

 src_file = sys.argv[1]
 dst_file = sys.argv[2]
 out_bands = 3

 def MakeColor(z_value):
  '''LUT for color ramp. Keys are pixel values, hash values are RGB
 triplets.'''
  color_dict = {
0:[102,0,255],
   25:[20,82,255],
   50:[0,194,224],
   75:[0,255,122],
  100:[41,255,0],
  125:[204,255,0],
  150:[245,194,0],
  175:[224,118,0],
  200:[168,46,0],
  225:[105,0,0],
  250:[64,0,0],
  255:[0,0,0]}

  key_list = color_dict.keys()
  key_list.sort()
  while len(key_list)  0:
last_val = key_list[-1]
if z_value = last_val:
  return color_dict[last_val]
else:
  key_list.remove(last_val)

 # Print some info
 print Creating %s % (dst_file)

 # Open source file
 src_ds = gdal.Open( src_file )
 src_band = src_ds.GetRasterBand(1)

 # create destination file
 dst_driver = gdal.GetDriverByName('GTiff')
 dst_ds = dst_driver.Create(dst_file, src_ds.RasterXSize,
 src_ds.RasterYSize, out_bands, 

[gdal-dev] Recommendation for what vector format to use?

[Roger André]

Hi All,

I'm about to implement a new workflow that uses gdal_rasterize and it dawned
on me that A) I wasn't limited to using shapefiles, and B) that it would be
a heck of alot more convenient to use a binary data format that stored all
the information in a single file.

Is there such a thing?  Do you all have any recomendations for what a good
interchange format is?

Thanks,

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

Re: [gdal-dev] need suggestion on Linux OS

[Roger André]

I have been running a fairly complete open source GIS stack on a couple of
Ubuntu 8.04 machines for about 6 months.  I have had no problems getting
anything installed from source on them, and use the supporting libs that are
available via apt-get almost exclusively.  If anyone wants to see the
program manifest and the install steps, I'd be happy to send it to them.

Roger
--

On Wed, Apr 22, 2009 at 9:27 AM, pcr...@pcreso.com wrote:


 Hi.

 I suggest you look at OpenSuse. That's my preferred FOSS GIS platform for a
 few years now.

 The guys that maintain the OpenSUuse GIS repository are really efficient 
 helpful. The couple of times I've wanted an updated package there I've
 emailed them  it's been done in a couple of days.

 I've also compiled some from CVS, when I needed the very latest  haven't
 had any problems.

 The repository (v11.1, 32 bit) is at:

 http://download.opensuse.org/repositories/Application:/Geo/openSUSE_11.1/i586/

 but it includes packages for other versions  other distros,
 http://download.opensuse.org/repositories/Application:/Geo/


 Works well  I prefer running Suse to Debian variants... But both work 
 it's nice to have the choice.


 Cheers,

   Brent wood


 --- On Thu, 4/23/09, John Donovan j.dono...@virtalis.com wrote:

  From: John Donovan j.dono...@virtalis.com
  Subject: RE: [gdal-dev] need suggestion on Linux OS
  To: Gong, Shawn (Contractor) shawn.g...@drdc-rddc.gc.ca, gdal-dev
 gdal-dev@lists.osgeo.org, openev-disc...@lists.sourceforge.net
  Date: Thursday, April 23, 2009, 2:41 AM
  I'm building the latest gdal
  (literally as I type!), on Kubuntu Hardy
  with no problems so far.
 
  -JD
 
   -Original Message-
   From: gdal-dev-boun...@lists.osgeo.org
 
   [mailto:gdal-dev-boun...@lists.osgeo.org]
  On Behalf Of Gong,
   Shawn (Contractor)
   Sent: 22 April 2009 15:24
   To: gdal-dev; openev-disc...@lists.sourceforge.net
   Subject: [gdal-dev] need suggestion on Linux OS
  
   hi list,
  
   I'd appreciate suggestions on Linux OS. We'd like to
  build
   newest gdal, as well as OpenEV2 using GTK2.
  
   Have tried CentOS. But it libraries seem too
  conservative.
  
   Is Fedora 9 a better choice?
  
   thanks,
  
   Shawn
  
  
  
  __
   This email has been scanned for Virtalis Ltd by the
   MessageLabs Email Security System.
  
  __
  
  
 
  __
  This email has been scanned for Virtalis Ltd by the
  MessageLabs Email Security System.
  __
  ___
  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] SRTM data and Pixel is area vs. Pixel is point

[Roger André]

Morning,

I have a bit of a conundrum which I'm trying to find an answer for,
and which is prompting a certain amount of contentious debate at work.
 I've downloaded a 1-deg x 1-deg SRTM tile from the GLCF site,
SRTM_f03_n007e000.tif.  Now from its name, this tile should have its
origin at N 7, E 0.

When I view the metadata for the tile, I see a helpful note that says,
Metadata: AREA_OR_POINT=Area

I assume from this that the data is using PixellsArea Raster Space,
as defined at http://www.remotesensing.org/geotiff/spec/geotiff2.5.html.
 This definition states that The PixelIsArea raster grid space R,
which is the default, uses coordinates I and J, with (0,0) denoting
the upper-left corner of the image.  So my guess would be that the UL
corner coordinates of this raster should be (0, 8).  However, that is
not the case.

 - When I run gdalinfo against the tile, I get the following info:
  Upper Left  (  -0.0004167,   8.0004167)

- When I create a worldfile for the tile with gdal_translate, I get
the following pixel coords reported in it:
  -0.00, 8.33

- Also, when I create a shapefile of this tile's extents using
gdaltindex, I get the following extents reported:
Extent: (-0.000417, 6.999583) - (1.000417, 8.000417)

So my question is this, are we really all working with SRTM data that
is offset by 1/2 a pixel because of a discrepancy between
point-vs-area pixel registration?  This is coming to a head for us
because we need to do some raster analysis, and we're having a really
hard time understandind why the data keeps spanning meridians and
parallels.

Would appreciate any clarification or insights.

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

Re: [gdal-dev] SRTM data and Pixel is area vs. Pixel is point

[Roger André]

Hi Frank,

I think the problem we're having is with this, It does imply a single
row/column of overlap with adjacent grid squares but I believe that is
typical of distributed SRTM data.  I guess we were expecting there to
be zero overlap and for the tiles to line up edge-to-edge.

Thanks,

Roger
--

On Fri, Apr 24, 2009 at 9:59 AM, Frank Warmerdam warmer...@pobox.com wrote:
 Roger André wrote:

 Morning,

 I have a bit of a conundrum which I'm trying to find an answer for,
 and which is prompting a certain amount of contentious debate at work.
  I've downloaded a 1-deg x 1-deg SRTM tile from the GLCF site,
 SRTM_f03_n007e000.tif.  Now from its name, this tile should have its
 origin at N 7, E 0.

 When I view the metadata for the tile, I see a helpful note that says,
 Metadata: AREA_OR_POINT=Area

 I assume from this that the data is using PixellsArea Raster Space,
 as defined at http://www.remotesensing.org/geotiff/spec/geotiff2.5.html.
  This definition states that The PixelIsArea raster grid space R,
 which is the default, uses coordinates I and J, with (0,0) denoting
 the upper-left corner of the image.  So my guess would be that the UL
 corner coordinates of this raster should be (0, 8).  However, that is
 not the case.

  - When I run gdalinfo against the tile, I get the following info:
  Upper Left  (  -0.0004167,   8.0004167)

 - When I create a worldfile for the tile with gdal_translate, I get
 the following pixel coords reported in it:
  -0.00, 8.33

 - Also, when I create a shapefile of this tile's extents using
 gdaltindex, I get the following extents reported:
 Extent: (-0.000417, 6.999583) - (1.000417, 8.000417)

 So my question is this, are we really all working with SRTM data that
 is offset by 1/2 a pixel because of a discrepancy between
 point-vs-area pixel registration?  This is coming to a head for us
 because we need to do some raster analysis, and we're having a really
 hard time understandind why the data keeps spanning meridians and
 parallels.

 Roger,

 I don't see a problem with the data.  It appears to have edge pixels
 that are properly centered on the lat/long grid points.  It does imply
 a single row/column of overlap with adjacent grid squares but I believe
 that is typical of distributed SRTM data.

 What, specifically, in the above seems wrong to you?

 Where we do often have contentious debates is when a raster is marked
 as PixelIsPoint.  In this case there is arguments about how it ought
 to be represented by GDAL and there are often complaints that the
 approach taken is off by half a pixel.  I think it would be much easier
 if everyone in the world just distributed data in pixelisarea mode - no
 risk of confusion and error.

 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] Help with first GDAL program

[Roger André]

Hi Leith,

This should help get you started.

Roger
--

On Sun, Apr 26, 2009 at 11:25 PM, Leith Bade le...@leithalweapon.geek.nzwrote:

 Hi,

 I am new to GDAL and would like some guidance on a simple task that I am
 writing a Python script for.

 Given a set of GeoTIFFs I would like to find the max/min lat/long extent of
 the set.

 This will be for calculating pixel dimensions of a theoretical
 gdal_merge.py type output, but for tiling into Google Earth.

 Thanks,
 Leith Bade
 le...@leithalweapon.geek.nz


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

#! /usr/bin/python

import sys
from osgeo import gdal

filename = sys.argv[1]

fh = gdal.Open( filename )
xsize = fh.RasterXSize
ysize = fh.RasterYSize
band_type = fh.GetRasterBand(1).DataType
projection = fh.GetProjection()
geotransform = fh.GetGeoTransform()

# randre - print geotransform info
print x pixel size =, geotransform[1]
print x rotation =, geotransform[2]
print y rotation =, geotransform[4]
print y pixel size =, geotransform[5]
print ul_x =, geotransform[0]
print ul_y =, geotransform[3]

ulx = geotransform[0]
lrx = ulx + geotransform[1] * xsize
uly = geotransform[3]
lry = uly + geotransform[5] * ysize




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

Re: [gdal-dev] free raster elevation model download

[Roger André]

http://srtm.csi.cgiar.org/
--

On Wed, May 6, 2009 at 4:30 PM, Matt Klaric m...@klaric.com wrote:



 On Wed, May 6, 2009 at 4:22 PM, mohwawang mohwaw...@yahoo.com wrote:


 Is there a website that provides free download for raster elevation model
 (DEM) that can be used by the gdal utility gdal_contour? Thanks.


 NASA:
 ftp://e0srp01u.ecs.nasa.gov/srtm/

 ___
 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] FOSS4G 2009 Registration Opens

[Roger André]

Just out of curiosity, how is the location of FOSS4G determined?

Thanks,

Roger
--

On Fri, May 8, 2009 at 6:22 AM, Frank Warmerdam warmer...@pobox.com wrote:

 Folks,

 FOSS4G is always the premier open source geospatial event of the
 year, and an important venue for MapServer and a variety of related
 projects.  I would encourage folks to attend, and to present on
 GDAL related activities.

 Frank

 ---


 2009 Free and Open Source Software for Geospatial (FOSS4G)
 October 20-23, Sydney, Australia

  http://2009.foss4g.org/

  Register prior to the July 31st Early Bird Deadline to
  save on your registration fees!

  Join 500 other geospatial professionals in
  exploring and explaining the world of
  OPEN SOURCE geospatial software.

  Workshops, labs, presentations, and networking,
  all in beautiful Sydney, Australia.

 SUBMIT A PRESENTATION

  Submit a presentation online at
  http://2009.foss4g.org/presentations/

  The deadline for submissions is June 1st, 2009.

  FOSS4G presentations are 25 minute talks, with 5 minute
  question and answer sessions at the end. Presentations
  cover the use or development of open source geospatial
  software. Anyone can submit a presentation proposal and
  take part in the conference as a presenter. More
  information is available on the presentations page on
  the website.

 EXHIBITOR AND SPONSORSHIP OPPORTUNITIES

  For information on exhibitor and sponsorship opportunities,
  see http://2009.foss4g.org/sponsorship/ or contact Cameron Shorter,
  Conference Chair via email at cameron.shor...@gmail.com.

  We hope to see you in Sydney, Australia in October!

  Yours,

  The FOSS4G 2009 Conference Committee


 ___
 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] FOSS4G 2009 Registration Opens

[Roger André]

Hi Jeff,

I have no other questions, and thanks for explaining how the selection
process works.

Roger
--

On Mon, May 11, 2009 at 11:13 AM, Jeff McKenna 
jmcke...@gatewaygeomatics.com wrote:

 Roger André wrote:

 Just out of curiosity, how is the location of FOSS4G determined?


 There is an OSGeo conference committee[1] which manages this decision. The
 general process is:

 - conference committee develops a Request for hosting Proposals document
 - RFP is brought to OSGeo Board for approval
 - RFP for hosting is released
 - the committee analyzes the proposals on the criteria listed in the RFP[2]
 - the committee's decision is then brought to the OSGeo Board for approval.

 If you have more questions please join the conference_dev mailing list[3].
  Thanks.

 -jeff


 [1] http://wiki.osgeo.org/wiki/Conference_Committee
 [2] http://www.osgeo.org/conference/rfp
 [3] http://lists.osgeo.org/mailman/listinfo/conference_dev




 --
 Jeff McKenna
 FOSS4G Consulting and Training Services
 http://www.gatewaygeomatics.com/


 ___
 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] Relax KML driver about xmlns

[Roger André]

Are you doing any version-dependent actions, or just checking for one of
three versions?  If not, then yes, I'd say take the namespace requirement
out.

Roger
--

On Tue, May 19, 2009 at 4:30 PM, Mateusz Loskot mate...@loskot.net wrote:

 Hi,

 Would it make sense to relax a bit the KML driver here [1] so it
 accepts kml root without xmlns attribute?

 One user reported to me that he is not able to read this file [2]
 using OGR KML driver and I found this is because of missing xmlns in
 kml root element. The KML was produced by some unknown
 software shipped with some GPS logger.

 I'm wondering if there are more KML generators in use that
 produce kml element without xmlns attribute.


 [1]

 http://trac.osgeo.org/gdal/browser/trunk/gdal/ogr/ogrsf_frmts/kml/kml.cpp#L268
 [2] http://mateusz.loskot.net/tmp/gdal/no-xmlns.kml

 Best regards,
 --
 Mateusz Loskot, http://mateusz.loskot.net
 Charter Member of OSGeo, http://osgeo.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] gdalwarp and Molleweide projections

[Roger André]

Hi Matthew,

It would be useful to see what your source image looks like.  I have an
image that covers essentially global extents, and which is in EPSG:4326.  I
did a reprojection using gdalwarp in GDAL 1.6.1:

$ gdalwarp -s_srs EPSG:4326 -t_srs +proj=moll +lon_0=0 infile outfile

The results that look reasonable, or at least I see no horizontal flipping
or other artifacts in the result.

Any chance you could share a reduced resolution version of the data with the
list?

Roger
--



On Sun, Jul 5, 2009 at 11:04 AM, Matthew Vavrek mattvav...@hotmail.comwrote:

 I have just started trying to use gdalwarp and gdal_translate to reproject
 some world maps. I have managed to achieve some limited success with some
 projections, but I seem to get strange results with my images, especially
 with Mollweide. The Molleweide gives a horizontally flipped image, sometimes
 with an unreversed strip in the middle, and the top and bottom of the image
 seem stretched rather than compressed as I would expect. The code I used to
 translate is:

 gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:54009 input.tif
 moll-output.tif

 I also tried with -t_srs '+proj=moll +lon_0=0'.
 My initial image was a plain jpg, and I have tried creating a world file
 and using gdal_translate to georeference the image. I also tried using:

 gdalwarp s_srs EPSG:4326 t_srs EPSG:4326 input.tif output.tif

 to see if that would help, which it oddly did for some of the other
 projections I was trying (like polyconic). I have even tried different world
 map images, yet still get these strange results. I can't figure out what I'm
 doing wrong to get these totally bizarre images, so any help would be
 appreciated.

 Thanks
 Matthew
 ___
 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] python gdal installation issues

[Roger André]

try adding the path to libgdal.so.1 in /etc/ld.so.conf

$ find /usr -name libgdal.so.1
/usr/local/lib/libgdal.so.1

$ sudo vi /etc/ld.so.conf
-- add
/usr/local/lib

$ ldconfig

$ /sbin/ldconfig -p | grep gdal
libgdal.so.1 (libc6) = /usr/local/lib/libgdal.so.1
libgdal.so (libc6) = /usr/local/lib/libgdal.so
--



On Wed, Jun 10, 2009 at 2:47 PM, Benjamin R Welton 
benjamin.r.wel...@wmich.edu wrote:

 Hey all,

Im having a bit of an issue using the python gdal api's. I keep getting
 an error in python saying that it is unable to find libgdal.so.1. I am not
 having any issues using the various gdal binaries however (gdalinfo works).
  i looked at PYPI and  changed the gdal-config location setting in the
 setup.cfg file which allowed me to install the python module but when i
 attempt to use it i get the libgdal.so.1 file missing  error.

   Anyone know of a way to fix this issue? I am compiling gdal from source
 on CentOS 5.1 using python 2.4.3.

 Thanks,
 Ben

 ___
 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] Writing georeferencing and projection info to NetCDF with Python Scientific.IO.NetCDF

[Roger André]

Hi all,

After a considerable delay due to work on other projects, I would like to
get back to exploring the capability of the GDAL netCDF driver.  To that
end, I am trying to create some netCDF files in which I can control variable
names, number of dimensions etc.  I have been exploring the
Scientific.IO.NetCDF
Python module and am able to create generic netCDF files, but they lack
georeferencing and projection info.

Can someone out there point me to some documentation, or cut-and-paste some
code, which shows how this info could be added?

Thanks in advance,

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

[gdal-dev] Status of TMS minidriver in gdal_wms

[Roger André]

Hi All,

I'd like to use GDAL to access tiles that have been generated by TileCache.
It appears this should be possible in GDAL 1.7 via the TMS service.  I
*think* the example at
http://www.gdal.org/frmt_wms_openstreetmap_tms.xmlshows this.

I'd appreciate it if someone could clarify whether my expectations are
reasonable regarding this, and whether this functionality can be enabled in
a pre-1.7 version of GDAL.

Thanks,

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

[gdal-dev] Re: Status of TMS minidriver in gdal_wms

[Roger André]

I was able to get this to work using the patch checked in from the Toronto
code sprint. However, I still have questions.  I would like to take
advantage of the fact that I already have pre-rendered tiles in my
TileCache, but it seems that MapServer is still having to render all of the
requests that I'm making through the WMS service.  How can I use my
pre-existing TileCache tiles as the data source?  Do I need to match the
projection that they are created in when I make the request?  Do I need to
set resolution levels in the GDAL_WMS file?  I followed the syntax in the
OSM GDAL_WMS sample to make the requests work.

Thanks,

Roger
--

On Wed, Aug 19, 2009 at 3:18 PM, Roger André ran...@gmail.com wrote:

 Hi All,

 I'd like to use GDAL to access tiles that have been generated by
 TileCache.  It appears this should be possible in GDAL 1.7 via the TMS
 service.  I *think* the example at
 http://www.gdal.org/frmt_wms_openstreetmap_tms.xml shows this.

 I'd appreciate it if someone could clarify whether my expectations are
 reasonable regarding this, and whether this functionality can be enabled in
 a pre-1.7 version of GDAL.

 Thanks,

 Roger
 --

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

Re: [gdal-dev] Shortening execution time of Python script

[Roger André]

So the problem is that ArcGIS Server is expecting the data to be ready
(unzipped), but due to the large size of the data set, the unzip has not
finished?

Roger
--

On Wed, Sep 9, 2009 at 3:54 AM, David Shi davidg...@yahoo.co.uk wrote:

 I have a Python script that automatically downloads zip files containing
 large datasets from another server and then unzips the files to further
 process the data.

 It has been used as a geoprocessor of ArcGIS Server.

 The script works fine when two datasets each has several kilobytes size,
 but the script stops half way when datasets were about 11,000KBytes.

 I think that the execution time is too long and ArcGIS Server just simply
 killed the process.

 What actions can I try to reduce the execution time?

 ArcGIS Server only works on the basis of 32 bits and I was told that the
 maximum memory it can utilise is 4 MBytes.

 I should be grateful if someone can make suggestions/recommendations.

 Sincerely,

 David


 ___
 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: gdal_rasterize with a NetCDF file?

[Roger André]

Hey Scott,

You could extract that layer from the netCDF, burn it with gdal_rasterize as
a geotiff, convert the geoTiff into a GMT grid, and then possibly add that
back into the netCDF (somehow).

Sounds like a pain though.

Roger
--

On Mon, Sep 21, 2009 at 11:49 AM, Scott Lewis scott.le...@nsidc.org wrote:

 Ah, thanks.  I must have missed that.  Looks like I'll have to find
 another way to accomplish what I'm trying to do then!

 Thanks for your help!

 Scott

 Hermann Peifer wrote:
  The driver doesn't support updating netCDF files. gdalinfo --formats
 says:
 
   GTiff (rw+): GeoTIFF
   (...)
   netCDF (rw): Network Common Data Format
 
  The + indicates update support, which is obviously missing for netCDF
  format.
 
  Hermann
 
 
   Original Message  
  Subject: Re:gdal_rasterize with a NetCDF file?
  From: Scott Lewis scott.le...@nsidc.org
  To: Date: 21/09/2009 20:37
 
  Evan,
 
  Thanks, that seems to be getting me in the right direction: I'm no
  longer getting the segfault error, or ERROR 5 at all.
 
  However, now it's giving me a different error.  The first line has the:
  0...10...20...30...40...50...60...70...80...90...100 - done.
  message, like I would expect from when I do this with GeoTIFF files.
  But following that I am getting a whole bunch of this error, repeated
  many times:
 
  ERROR 6: WriteBlock() not supported for this dataset.
 
 
  Doing a search online, I'm getting the impression that this might mean
  that gdal_rasterize doesn't support writing to a NetCDF file?  Is this
  correct, or does the problem lie elsewhere with my calling of the
  utility.
 
  Here's the command line I'm using now:
 
  /usr/local/FWTools-2.0.6/bin_safe/gdal_rasterize -i -burn 0 -l
  OGRGeoJSON dummy.json
  'NETCDF:Albedo.nc:Albedo_with_1400_Local_Time_of_Measurement'
 
 
  Thanks again for the help!
 
  Scott
 
 
  Even Rouault wrote:
  Selon Scott Lewis scott.le...@nsidc.org:
 
  Scott,
 
  Anytime you want to access to a subdataset with a gdal utility, you
  must use the
  value of the relevant SUBDATASET_xxx_NAME metadata item, in that
  instance
  'NETCDF:Albedo.nc:Albedo_with_1400_Local_Time_of_Measurement'
 
  The segmentation fault however wasn't expected, so I commited fixes
  per ticket
  http://trac.osgeo.org/gdal/ticket/3146 to error out more nicely.
 
  Best regards
 
  Hi,
 
  I'm still fairly new at using GDAL, but I am having trouble burning a
  vector image onto a NetCDF file.  I've had success with using
  gdal_rasterize to burn a vector onto a GeoTIFF file, but when I try to
  burn the same vector onto a NetCDF file, I get an error.
 
  Based on the error, I'm suspecting it's because it doesn't know how to
  find the right band in the NetCDF File (I want the
  Albedo_with_1400_Local_Time_of_Measurement subdataset), but although
  I've looked at the documentation, I'm not sure how to specify this to
  gdal_rasterize.
 
  Any help would be much appreciated!  Below is the info.
 
 
 
  Here is the command I am trying to run:
 
  /usr/local/FWTools-2.0.6/bin_safe/gdal_rasterize -i -burn 0 -l
  OGRGeoJSON dummy.json Albedo.nc
 
 
  And the error I get:
 
  ERROR 5: GDALDataset::GetRasterBand(1) - Illegal band #
 
  /usr/local/FWTools-2.0.6/bin_safe/gdal_rasterize: line 9: 26950
  Segmentation fault  $FWTOOLS_HOME/bin/`basename $TARGET` $@
 
 
 
  Here is the gdalinfo for Albedo.nc:
 
  Driver: netCDF/Network Common Data Format
  Files: Albedo.nc
  Size is 512, 512
  Coordinate System is `'
  Metadata:
NC_GLOBAL#Conventions=CF-1.4
NC_GLOBAL#institution=National Snow  Ice Data Center, Boulder, CO
NC_GLOBAL#title=AVHRR Polar Pathfinder Twice-Daily 5 km EASE-Grid
  Composites Albedo
NC_GLOBAL#source=See the title and references for this information
NC_GLOBAL#comment=Not set at this time
NC_GLOBAL#references=Documentation available at:
  http://nsidc.org/data/docs/daac/nsidc0066_avhrr_5km.gd.html
NC_GLOBAL#history=Mon, 21 Sep 2009 09:31:24: File created.
  Subdatasets:
 
 
 
 SUBDATASET_1_NAME=NETCDF:Albedo.nc:Albedo_with_1400_Local_Time_of_Measurement
 
SUBDATASET_1_DESC=[1x244x242]
  Albedo_with_1400_Local_Time_of_Measurement (8-bit integer)
SUBDATASET_2_NAME=NETCDF:Albedo.nc:latitude
SUBDATASET_2_DESC=[244x242] latitude (64-bit floating-point)
SUBDATASET_3_NAME=NETCDF:Albedo.nc:longitude
SUBDATASET_3_DESC=[244x242] longitude (64-bit floating-point)
  Corner Coordinates:
  Upper Left  (0.0,0.0)
  Lower Left  (0.0,  512.0)
  Upper Right (  512.0,0.0)
  Lower Right (  512.0,  512.0)
  Center  (  256.0,  256.0)
 
 
  And for completeness, the dummy.json file (it's basically a triangle):
 
  {
  type: FeatureCollection,
  features: [{
  type: Feature,
  properties: { BASIN_ID: 1.00, AREA_SQKM:
  1739.016470 },
  geometry: {
  type: Polygon,
  coordinates: [ [
  [ -200, -250 ],
  [ -250, -140 ],

Re: [gdal-dev] Python bindings

[Roger André]

Had a similar problem in the past that I resolved like this:

$ gdalinfo --formats
gdalinfo: error while loading shared libraries: libgdal.so.1: cannot
open shared object file: No such file or directory

-- Is the lib acessible outside of gdal?
$ /sbin/ldconfig -p | grep gdal

--nuffin...  Let's try adding the path to libgdal.so.1 in /etc/ld.so.conf

$ find /usr -name libgdal.so.1
/usr/local/lib/libgdal.so.1

$ sudo vi /etc/ld.so.conf
-- add
/usr/local/lib

$ ldconfig

$ /sbin/ldconfig -p | grep gdal
libgdal.so.1 (libc6) = /usr/local/lib/libgdal.so.1
libgdal.so (libc6) = /usr/local/lib/libgdal.so
---

On Sat, Oct 10, 2009 at 5:44 AM, Chaitanya kumar CH
chaitanya...@gmail.com wrote:
 Hi,

 I am having trouble running GDAL's Python bindings. I compiled the trunk
 with the option --with-python. Below is some output from my command line. I
 am running this on a newly installed system. It has only 2.6 version of
 Python and python-setuptools.

 --
 ubu...@ubuntu:~/gdal/2270/trunk/gdal$ export PYTHONPATH=/usr/local/
 ubu...@ubuntu:~/gdal/2270/trunk/gdal$ python
 Python 2.6.2 (release26-maint, Apr 19 2009, 01:58:18)
 [GCC 4.3.3] on linux2
 Type help, copyright, credits or license for more information.
 from osgeo import gdal
 Traceback (most recent call last):
  File stdin, line 1, in module
  File
 /usr/local/lib/python2.6/dist-packages/GDAL-1.6.0-py2.6-linux-x86_64.egg/osgeo/__init__.py,
 line 4, in module
    import _gdal
 ImportError: libgdal.so.1: cannot open shared object file: No such file or
 directory
 import gdal
 Traceback (most recent call last):
  File stdin, line 1, in module
  File
 /usr/local/lib/python2.6/dist-packages/GDAL-1.6.0-py2.6-linux-x86_64.egg/gdal.py,
 line 2, in module
    from osgeo.gdal import deprecation_warn
  File
 /usr/local/lib/python2.6/dist-packages/GDAL-1.6.0-py2.6-linux-x86_64.egg/osgeo/__init__.py,
 line 4, in module
    import _gdal
 ImportError: libgdal.so.1: cannot open shared object file: No such file or
 directory

 ubu...@ubuntu:~/gdal/2270/trunk/gdal$ ls /usr/local/lib/
 libgdal.a  libgdal.la  libgdal.so  libgdal.so.1  libgdal.so.1.13.0
 python2.6
 ubu...@ubuntu:~/gdal/2270/trunk/gdal$
 --

 Thank you.
 --
 Chaitanya kumar CH.


 ___
 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] open source tool

[Roger André]

This could be done with GRASS, I think.  First extract the LAS file
into an ASCII x,y,z file.  Then load into GRASS with r.in.xyz. Then
use one of the v.delaunay or v.voronoi GRASS modules.

Roger
--

2009/8/27 Riki Tiki mongoose_r...@mail.ru:
 hi,

 do you know any tools available for the conversion of LAS (LIDAR data) into 
 Voronoi diagrams or TIN (triangular irregular networks)? thanks a lot!

 or anyone can suggest how to do it in arcgis for example?

 thanks


 ___
 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] Setting nodata tag with Python API

[Roger André]

Hi All,

I'm using the Python API to build some GeoTIFFs and would like to set
the nodata value to - at creation time.  How do I do that?  Below
is a sample of my code

  def makeRaster(self):
out_drv = gdal.GetDriverByName('GTiff')
out_ds = out_drv.Create(self.raster_name, 4320, 2160 , 1, gdal.GDT_Float32)
out_ds.SetGeoTransform([-180, 0.080, 0, 90, 0, -0.080])
out_band = numpy.zeros([out_ds.RasterYSize, out_ds.RasterXSize])
for iY in range(out_ds.RasterYSize):
  for iX in range(out_ds.RasterXSize):
alloc_key = self.makeAllocKey(iY + 1, iX + 1)
pixel_value = self.checkPixel(alloc_key)
out_band[iY][iX] = pixel_value
out_ds.GetRasterBand(1).WriteArray(out_band)
out_ds = None

Thanks,

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

Re: [gdal-dev] trying to change regionator image resampling to nearest-neighbor

[Roger André]

Hi Mary Jo,

Sounds like you have started down a slightly different path to resolve
this by using gdal2tiles.py, but I'm still a bit curious about the
results you were seeing with the regionator.  Would you mind sending
me a couple screenshots that show the before and after images, with
haloing?  Your description makes me think that it is a jpg compression
problem, more than a resampling artifact.  It might be possible to
simply change the output format to PNG in the code, and eliminate the
problem.

Regards,

Roger
--



On Thu, Nov 12, 2009 at 11:32 AM, Mary Jo Brodzik brod...@nsidc.org wrote:

 Dear gdal-dev list members,

 I'm using a python tool called the regionator
 (http://code.google.com/p/regionator/wiki/ Welcome) on a large (14352x7440)
 png image.  The regionator uses the GDAL library to resample subsets of the
 original image into much smaller chunk images.  These images are set up at
 increasing levels of detail in a kml file that make the download of this
 information into a googleEarth client much more efficient.

 The png is color and includes transparency; its approximate resolution is
 500 meters. I basically have 4 colors in the original image, black, white,
 blue and yellow. The output images that regionator is producing are great,
 but when I zoom in to distances where the original pixel sizes are visible,
 it looks like the resampling did some averaging across the 4 colors, so I'm
 getting blue sections surrounded by a blue halo that bleeds into white, and
 yellow and blue values bleeding into each other, etc.

 I'd like to control the resampling, and force it to do nearest-neighbor
 rather than averaging to eliminate the halos.  It looks to me like the
 regionator is doing the resampling with calls to ReadRaster and WriteRaster
 in in a class called extractor.py, but I don't know enough about the GDAL
 libraries to know how to change them. Is it possible to force the resampling
 to be nearest neighbor?  extractor.py isn't terribly large or complicated,
 so I'm including it here (below).

 I've already tried the googleearth forum and the kml-developers forum, no
 help from those communities.  Python's new to me, and I've only used the
 GDAL command-line utilities up to now.

 Thanks in advance for any tips.

 Mary Jo

 Here's the regionator's extractor.py:

 
 Copyright (C) 2006 Google Inc.

 Licensed under the Apache License, Version 2.0 (the License);
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an AS IS BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 

 
 $URL: http://regionator.googlecode.com/svn/trunk/kml/extractor.py $
 $Revision: 200 $
 $Date: 2007-01-12 09:12:48 -0700 (Fri, 12 Jan 2007) $
 

  class Extractor

 Phase 3 of SuperOverlay creation.

 Extract tiles from an image

 Front-ends GDAL.  Input file and output driver must be
 some format known to GDAL.

 

 import gdal
 import tempfile
 import os


 class Extractor:

   class Extractor

  Create for Extractor the image file.
  Each tile is resampled to the specified size and written
  to a file of the given format.

  

  def __init__(self,imgfile,twid,tht,fmt,verbose=False):

    

    Args:
      imgfile: a gtiff
      twid,tht: extracted tile resample pixel width/height
      fmt: GDAL output driver name ('PNG','JPEG')
    

    self.__in_ds = gdal.Open(imgfile)
    self.__twid = twid
    self.__tht = tht
    self.__fmt = fmt
    self.__verbose = verbose

    self.__bands = self.__in_ds.RasterCount

    if self.__verbose:
      print 'Extractor %s %d bands' % (imgfile,self.__bands)

    # Intermediate work must be in GTiff (?)
    self.__gtiff_driver = gdal.GetDriverByName('GTiff')

    # Output driver/format is whatever the user specifies
    # XXX handle bad/wrong fmt's less gracelessly
    self.__o_driver = gdal.GetDriverByName(fmt)

  def Extract(self,x,y,wid,ht,basename):

     Extract a tile into a file

    The given tile is extracted, resampled and saved
    according to the tile pixel dimensions and format
    specified at object __init__().

    Args:
      x,y: pixel offset
      wid,ht: pixel dimensions

    Returns:
      True: complete success
      False: any failure

    

    # Get the tile's pixels resampled
    twid = self.__twid
    tht = self.__tht
    i_data = self.__in_ds.ReadRaster(x,y,wid,ht,buf_xsize=twid,buf_ysize=tht)

    # Have to Create out to GTiff first (?)
    (fd, tmpfile) = tempfile.mkstemp(suffix='GTiff')
    os.close(fd)
    o_ds = self.__gtiff_driver.Create(tmpfile,twid,tht,bands=self.__bands)
    o_ds.WriteRaster(0,0,twid,tht,i_data)

    # Save off using the specified driver
    filename = '%s.%s' % (basename, 

[gdal-dev] Can GDAL open binary buffer contents as input for MEM format?

[Roger André]

Hi All,

I'd like to know if I can take the results of a MapServer Python mapscript
'mapImage.getBytes()' operation, and feed it into GDAL as a data source for
an in-memory raster?  The MapServer docs state that the results of
'getBytes' in Python is a string of binary data.  I'm not sure if GDAL can
open this sort of object though.

For background purposes, the reason I would like to do this is so that I can
try replacing PIL with GDAL in TileCache for the cutting of metatiles.  PIL
isn't properly saving the image formats generated by MapServer, whereas GDAL
replicates them perfectly.

Thanks,

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

Re: [gdal-dev] error about gdal shared libraries from within GRASS

[Roger André]

Hi Isaac,

I see that you've updated LD_LIBRARY_PATH already, but perhaps you could try
adding the lib location to ls.so.conf as well.  See below for how I dealt
with similar error in the past.

--
$ gdalinfo --formats
gdalinfo: error while loading shared libraries: libgdal.so.1: cannot open
shared object file: No such file or directory

-- ruh roh
$ /sbin/ldconfig -p | grep gdal

--nuffin...  Let's try adding the path to libgdal.so.1 in /etc/ld.so.conf

$ find /usr -name libgdal.so.1
/usr/local/lib/libgdal.so.1

$ sudo vi /etc/ld.so.conf
-- add
/usr/local/lib

$ ldconfig

$ /sbin/ldconfig -p | grep gdal
libgdal.so.1 (libc6) = /usr/local/lib/libgdal.so.1
libgdal.so (libc6) = /usr/local/lib/libgdal.so

-- alright, so what about gdalinfo?
$ gdalinfo
Usage: gdalinfo [--help-general] [-mm] [-stats] [-nogcp] [-nomd]
[-noct] [-checksum] [-mdd domain]* datasetname


Best of luck,

Roger
--


On Fri, Jan 15, 2010 at 10:13 AM, Isaac Ullah isaac.ul...@asu.edu wrote:

 Hi list,

 I have recently been trying to compile and install 64 bit GRASS7 with
 64 bit GDAL1.6.3 on 64 bit Ubuntu 9.10. I have successfully compiled both
 packages (with no errors), but am getting this error upon GRASS startup when
 the GRASS gui loads and first tries to access GDAL:
 g.list: error while loading shared libraries: libgdal.so.1: cannot open
 shared object file: No such file or directory

 I have that library installed in a standard location:
  whereis libgdal.so.1
 libgdal.so: /usr/lib/libgdal.so.1 /usr/lib64/libgdal.so.1
 /usr/local/lib/libgdal.so /usr/local/lib/libgdal.so.1

 and even have made gdal.conf:
  cat /etc/ld.so.conf.d/gdal.conf
 /usr/local/lib

 And also, this file IS 64 bit:
  file /usr/local/lib/libgdal.so.1.13.3
 /usr/local/lib/libgdal.so.1.13.3: ELF 64-bit LSB shared object, x86-64,
 version 1 (SYSV), dynamically linked, not stripped

 Even after explicitly setting
 LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib before grass startup,
 libgdal.so.1 is not found, even though echo $LD_LIBRARY_PATH produces:
 /usr/local/grass7.0.svn-i686-pc-linux-gnu-28_11_2009/lib:/usr/local/lib

 Result of ldd /usr/local/grass7.0.svn-28_11_2009/bin/g.list shows the only
 problem: libgdal.so.1 = not found.

 Here are some additional ldd responses:
 ldd `which gdalinfo`
 linux-vdso.so.1 =  (0x7fff207ff000)
 libgdal.so.1 = /usr/local/lib/libgdal.so.1 (0x7fc141b91000)
 libgeos_c.so.1 = /usr/local/lib/libgeos_c.so.1 (0x7fc141975000)
 libgeos-3.2.0.so = /usr/local/lib/libgeos-3.2.0.so(0x7fc1415fa000)
  libsqlite3.so.0 = /usr/lib/libsqlite3.so.0 (0x7fc141372000)
 libodbc.so.1 = /usr/lib/libodbc.so.1 (0x7fc141112000)
 libodbcinst.so.1 = /usr/lib/libodbcinst.so.1 (0x7fc140f05000)
 libltdl.so.7 = /usr/lib/libltdl.so.7 (0x7fc140cfb000)
 libexpat.so.1 = /lib/libexpat.so.1 (0x7fc140ad2000)
 libxerces-c.so.28 = /usr/lib/libxerces-c.so.28 (0x7fc1404ff000)
 libjasper.so.1 = /usr/lib/libjasper.so.1 (0x7fc1402a8000)
  libjpeg.so.62 = /usr/lib/libjpeg.so.62 (0x7fc140083000)
 libhdf5-1.6.6.so.0 = /usr/lib/libhdf5-1.6.6.so.0 (0x7fc13fd6)
 libz.so.1 = /lib/libz.so.1 (0x7fc13fb49000)
 libmfhdf.so.4 = /usr/lib/libmfhdf.so.4 (0x7fc13f921000)
  libdf.so.4 = /usr/lib/libdf.so.4 (0x7fc13f66f000)
 libnetcdf.so.4 = /usr/lib/libnetcdf.so.4 (0x7fc13f43c000)
 libpq.so.5 = /usr/lib/libpq.so.5 (0x7fc13f214000)
 libpthread.so.0 = /lib/libpthread.so.0 (0x7fc13eff8000)
  librt.so.1 = /lib/librt.so.1 (0x7fc13edf)
 libdl.so.2 = /lib/libdl.so.2 (0x7fc13ebec000)
 libcurl-gnutls.so.4 = /usr/lib/libcurl-gnutls.so.4
 (0x7fc13e9aa000)
 libstdc++.so.6 = /usr/lib/libstdc++.so.6 (0x7fc13e69a000)
  libm.so.6 = /lib/libm.so.6 (0x7fc13e416000)
 libgcc_s.so.1 = /lib/libgcc_s.so.1 (0x7fc13e1ff000)
 libc.so.6 = /lib/libc.so.6 (0x7fc13de9)
 libicuuc.so.40 = /usr/lib/libicuuc.so.40 (0x7fc13db48000)
  libicudata.so.40 = /usr/lib/libicudata.so.40 (0x7fc13cc03000)
 libssl.so.0.9.8 = /lib/libssl.so.0.9.8 (0x7fc13c9b5000)
 libcrypto.so.0.9.8 = /lib/libcrypto.so.0.9.8 (0x7fc13c62e000)
 libkrb5.so.3 = /usr/lib/libkrb5.so.3 (0x7fc13c376000)
 libcom_err.so.2 = /lib/libcom_err.so.2 (0x7fc13c172000)
 libgssapi_krb5.so.2 = /usr/lib/libgssapi_krb5.so.2
 (0x7fc13bf44000)
 libcrypt.so.1 = /lib/libcrypt.so.1 (0x7fc13bd0b000)
 libldap_r-2.4.so.2 = /usr/lib/libldap_r-2.4.so.2 (0x7fc13bac2000)
 /lib64/ld-linux-x86-64.so.2 (0x7fc142692000)
 libidn.so.11 = /usr/lib/libidn.so.11 (0x7fc13b88f000)
  liblber-2.4.so.2 = /usr/lib/liblber-2.4.so.2 (0x7fc13b681000)
 libgnutls.so.26 = /usr/lib/libgnutls.so.26 (0x7fc13b3df000)
 libgcrypt.so.11 = /lib/libgcrypt.so.11 (0x7fc13b167000)
 libk5crypto.so.3 = /usr/lib/libk5crypto.so.3 

Re: [gdal-dev] Can GDAL open binary buffer contents as input for MEM format?

[Roger André]

Hi Frank,

Thanks for the example.  How can I specify alpha in that, for one of the
quantized 8-bit PNG's (with alpha) that MapServer produces?

They look like this:

Metadata:
  Software=UMN Mapserver
Corner Coordinates:
Upper Left  (0.0,0.0)
Lower Left  (0.0,  256.0)
Upper Right (  256.0,0.0)
Lower Right (  256.0,  256.0)
Center  (  128.0,  128.0)
Band 1 Block=256x1 Type=Byte, ColorInterp=Palette
  NoData Value=0
  Color Table (RGB with 85 entries)
0: 0,0,0,0
1: 251,0,0,64
2: 250,0,0,54
3: 253,0,0,114
4: 254,0,0,176
5: 254,0,0,224
6: 254,0,0,228
7: 254,0,0,252
8: 253,0,0,168
9: 254,0,0,170
   10: 254,0,0,234
   11: 254,0,0,246
   12: 254,0,0,230
   13: 252,0,0,194
   14: 253,0,0,232
   15: 254,0,0,248
   16: 254,0,0,226
   17: 253,0,0,252
   18: 254,0,0,244
   19: 254,0,0,240
   20: 253,0,0,248
   21: 254,0,0,192
   22: 250,0,0,98
   23: 252,0,0,98
   24: 253,0,0,236
   25: 251,0,0,144
   26: 253,0,0,208
   27: 252,0,0,252
   28: 253,0,0,250
   29: 253,0,0,226
   30: 254,0,0,212
   31: 253,0,0,242
   32: 250,0,0,112
   33: 253,0,0,238
   34: 253,0,0,222
   35: 254,0,0,208
   36: 253,0,0,234
   37: 254,0,0,236
   38: 254,0,0,232
   39: 252,0,0,150
   40: 253,0,0,240
   41: 251,0,0,114
   42: 253,0,0,244
   43: 254,0,0,222
   44: 231,0,0,255
   45: 202,0,0,255
   46: 249,0,0,255
   47: 96,0,0,255
   48: 63,0,0,255
   49: 229,0,0,255
   50: 89,0,0,255
   51: 248,0,0,255
   52: 192,0,0,255
   53: 230,0,0,255
   54: 88,0,0,255
   55: 176,0,0,255
   56: 225,0,0,255
   57: 206,0,0,255
   58: 64,0,0,255
   59: 53,0,0,255
   60: 113,0,0,255
   61: 169,0,0,255
   62: 183,0,0,255
   63: 252,0,0,255
   64: 184,0,0,255
   65: 251,0,0,255
   66: 247,0,0,255
   67: 189,0,0,255
   68: 224,0,0,255
   69: 250,0,0,255
   70: 244,0,0,255
   71: 243,0,0,255
   72: 227,0,0,255
   73: 241,0,0,255
   74: 253,0,0,255
   75: 246,0,0,255
   76: 167,0,0,255
   77: 168,0,0,255
   78: 254,0,0,255
   79: 242,0,0,255
   80: 245,0,0,255
   81: 44,0,0,255
   82: 161,0,0,255
   83: 0,0,0,255
   84: 255,0,0,255


Roger
--

On Sat, Jan 16, 2010 at 1:26 PM, Frank Warmerdam warmer...@pobox.comwrote:

 Roger André wrote:

 Hi All,

 I'd like to know if I can take the results of a MapServer Python mapscript
 'mapImage.getBytes()' operation, and feed it into GDAL as a data source for
 an in-memory raster?  The MapServer docs state that the results of
 'getBytes' in Python is a string of binary data.  I'm not sure if GDAL can
 open this sort of object though.

 For background purposes, the reason I would like to do this is so that I
 can try replacing PIL with GDAL in TileCache for the cutting of metatiles.
  PIL isn't properly saving the image formats generated by MapServer, whereas
 GDAL replicates them perfectly.


 Roger,

 I am not aware of a way to utilize a raw buffer as a mem dataset from
 Python (though this is doable from C++).  However, you can just write
 the buffer to a mem dataset.  This code snippet from the autotest suite
 shows something like this:

###
# Setup dataset
drv = gdal.GetDriverByName('MEM')
gdaltest.mem_ds = drv.Create( 'mem_1.mem', 50, 3, gdal.GDT_Int32, 1 )
ds = gdaltest.mem_ds

raw_data = array.array('f',list(range(150))).tostring()
ds.WriteRaster( 0, 0, 50, 3, raw_data,
buf_type = gdal.GDT_Float32,
band_list = [1] )

 It is slightly more tricky with what I assume is a pixel interleaved 3
 band image you would have.  Without having tested it, something like this
 might work:

  mem_ds = gdal.GetDriverByName('MEM').Create('x',200,300,gdal.GDT_Byte,3)
  mem_ds.WriteRaster( 0, 0, 200, 300,
  mapImage.getBytes(), 200, 300, gdal.GDT_Byte,
  band_list=[1,2,3],
  buf_pixel_space=3,
  buf_line_space=3*200,
  buf_band_space=1)

 Good luck,


 --

 ---+--
 I set the clouds in motion - turn up   | Frank Warmerdam,
 warmer...@pobox.com
 light and sound - activate the windows | 
 http://pobox.com/~warmerdamhttp://pobox.com/%7Ewarmerdam
 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] Can GDAL open binary buffer contents as input for MEM format?

[Roger André]

Hi Frank,

Ok, that makes sense.  So for this to work, I will have to know (somehow)
what type of image mapserver is sending, and then instantiate a MEM dataset
with the correct parameters.  Problem with that is that I can't think of a
good way to read what those parameters are out of the binary string.  It's
no problem when I know exactly what data is coming in, as is the case for my
very specific need, but it means that swapping in GDAL for PIL as a general
solution might not be possible.

I'll play around with the examples you've sent, and will investigate whether
there are ways to read the string and determine what data type it is.

One last question. How does a utility like gdal_translate figure out what
the input data type is?

Thanks again,

Roger
--


On Sat, Jan 16, 2010 at 2:34 PM, Frank Warmerdam warmer...@pobox.comwrote:

 Roger André wrote:

 Hi Frank,

 Thanks for the example.  How can I specify alpha in that, for one of the
 quantized 8-bit PNG's (with alpha) that MapServer produces?


 Roger,

 In a way this is simplier.  It is just a one band GDT_Byte file, so more
 like the first example I gave.  The palette can be handled separately,
 and assigned to the MEM dataset a bit like:

gdaltest.test_ct_data = [ (255,0,0), (0,255,0), (0,0,255),
 (255,255,255,0)]

gdaltest.test_ct = gdal.ColorTable()
for i in range(len(gdaltest.test_ct_data)):
gdaltest.test_ct.SetColorEntry( i, gdaltest.test_ct_data[i] )
mem_ds.GetRasterBand(1).SetColorTable( gdaltest.test_ct )

 Best regards,

 --

 ---+--
 I set the clouds in motion - turn up   | Frank Warmerdam,
 warmer...@pobox.com
 light and sound - activate the windows | 
 http://pobox.com/~warmerdamhttp://pobox.com/%7Ewarmerdam
 and watch the world go round - Rush| Geospatial Programmer for Rent


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

[gdal-dev] Using RasterIO in Python to crop images

[Roger André]

Hi all,

I'd like to do some cropping based on pixel coordinates in a script I'm
writing.  I'm familiar with using the srcwin option to gdal_translate, but
have not done this using the API before.  Reading the API docs and some of
the previous postings to this list, it appears that the RasterIO function is
exposed in the Python API via band.ReadRaster().  I'd like some confirmation
that I'm understanding the docs correctly.

Assuming that I have an 8bit, paletted image named foo.png, if I wanted to
extract a 50 x 50 pixel tile that had it's origin at pixel coords (0,0), I
would do something like this:

in_ds = gdal.Open('foo.png')
band = dataset.GetRasterBand(1)
tile_data = band.ReadRaster( 0, 0, 50, 50, 2500, 1, GDT_Byte )

and then I would imagine I have to create the output dataset and write
tile_data to it.

Is there a way to do specify the output dataset in such a way that it
matches all the params of the input data, but just contains the smaller data
area?

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

[gdal-dev] How to calculate total change in headings?

[Roger André]

Hi all,

Not a GDAL question at all, but one which might amuse some of you, and which
has been unsolveable in a satisfactory fashion for me.

Given a list of compass headings in degrees, how can you reliably calculate
the total change in direction?

Here's an example:  [355.5, 10.0, 5.1]

This seemed like a simple thing to figure out, but seems to be more
complicated than I had originally thought.

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

Re: [gdal-dev] Raster data extraction

[Roger André]

Hi Brent,

You can write a program that takes in the coordinates of your sample
points as input, and then reads the pixel values in each band of the
image at that location.

Not sure if that's what you're looking for, but it works.

Roger
--

On Mon, Apr 19, 2010 at 2:36 PM,  pcr...@pcreso.com wrote:
 Hi,

 Is there a simple way (using gdal) of retrieving the colour spectra at 
 specified sample points in georeferenced images?

 I figure GRASS has some tools to do this, as does GMT, but was wondering if 
 GDAL could support this functionality directly.

 Lots of points/images to work through.

 Thanks,

  Brent Wood
 ___
 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] SRTM/hgt to TIFF with alpha

[Roger André]

Hi David,

Can you clip out a small section of your DEM that contains both ocean and
land and send it to me?  Your question sounds interesting, and I'd like to
take a look at the data before I comment on it.

Thanks,

Roger
--

On Wed, Apr 7, 2010 at 9:08 PM, David dbainbri...@gmail.com wrote:

 Hello.  I  am having a terrible time trying to figure out how to convert
 what should be value/height of 0 in the SRTM/hgt data file to the alpha
 layer in a tiff.  I was experimenting with making ocean areas transparent.
  I tried the following with no success:

 gdalwarp -of GTiff -co TILED=YES  -dstalpha -srcnodata 0 -wt Float32 -ot
 Float32 -wo SAMPLE_STEPS=100 -multi N40W125.hgt test.tiff

 gdalinfo -mm gives:
 Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
Computed Min/Max=-32768.000,1233.000
 Band 2 Block=256x256 Type=Float32, ColorInterp=Alpha
Computed Min/Max=0.000,255.000

 Since this may affect the hillshading for any land areas inland I was also
 going to attempt the same affect by masking the ocean after the hillshading.

 The resulting tiff in either case would be overlaid onto another layer that
 would have any land areas filled in with the resulting color that would have
 been used with height 0 in the hgt file.  I am doing this because I have a
 coastline shapefile that represents land areas but not ocean areas.

 Is my approach a good idea and/or what am I missing to mask the ocean areas
 into the alpha layer?
 Thanks for any help.

 ___
 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] create a png file

[Roger André]

Have you looked at http://www.gdal.org/gdal_tutorial.html?

Attached is a Python code snippet that shows how to assign values to pixels
in an image. To do what you're asking, you need to create the output raster,
and then write out the values in your array into the correct pixel position.

I assume that you need to modify the values in the DEM, and not just simply
convert the image from one format to another?

Roger
--

2010/4/5 weixj2003ld weixj200...@163.com

 Now,I can read data from a DEM file into an array,but I do not know how to
 use it to create a png file?
 Thk u.
 BTW,Some one tell me that freeimage can do it,but  I do not know.




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

#! /usr/bin/env python

Create an image with gdal using hard-coded georeferencing and extents.

 tile_span = 5_deg x 5_deg
x_span = 300 pixels
y_span = 300 pixels
pixel_size = .0166
ul_x = -60
ul_y = -5


import gdal
import Numeric

# .GetDriverByName(gdal_format)
# .create(filename, x_span, y_span, bands, data_type)
# .SetGeoTransform([ul_x, pixel_x, rotation_x, ul_y, rotation_y, pixel_y]) 

out_drv = gdal.GetDriverByName('GTiff')
out_ds = out_drv.Create('step1.tif', 300, 300 , 1, gdal.GDT_Byte)
out_ds.SetGeoTransform([-60, .0166, 0, -5, 0, -.0166])
out_band = Numeric.zeros([out_ds.RasterYSize, out_ds.RasterXSize])

# walk through the arrays and assign a value for each pixel
for iY in range(out_ds.RasterYSize):
  for iX in range(out_ds.RasterXSize):
out_band[iY][iX] = 0
  
out_ds.GetRasterBand(1).WriteArray(out_band)

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

Re: [gdal-dev] NetCDF geotransform - pixel centre/upper-left issue

[Roger André]

From my experience working with NetCDF files in GDAL, your best option is to
write your own tools to work between NetCDF and other raster formats.
Specifically, I would use the API to read the NetCDF, and explicitly define
the georeferencing of your output file using whatever logic you know to be
correct for your input data.  All of the variables in the NetCDF file are
accessible via GDAL, but in my experience the built-in utilities do not
always use them in the way you want.  I too believe this is a failure of the
NetCDF format, and not in GDAL.

Best of luck,

Roger
--

On Mon, Apr 5, 2010 at 11:26 PM, Michael Sumner mdsum...@gmail.com wrote:

 I believe that NetCDF generally uses cell centres to specify pixel
 coordinates, and there is no way to specify cell corners except by
 adding metadata about the actual dimensions boundaries (this is more
 problematic again if coordinates are not regular since you cannot
 assume a regular half cell).

 My understanding is that GDAL assumes cell centres for NetCDF, and
 from the data I've seen this is reasonable. The data model in NetCDF
 is more general than GDAL in some ways and it's simply not possible to
 answer this issue completely.  The unhappy news is that there is no
 easy answer but IMO GDAL is not incorrect here.

 In that example file the first longitude is 0 and all subsequent
 values are regularly spaced from there up to one cell less than 360,
 so the cell bounds at [0,360] are obvious - but this is not always the
 case. The latitudes are not regular and what was intended cannot
 always be determined in many cases (although in this file the
 boundaries are fairly obviously at [-90,90] offset from the min/max
 values of -87.8638  87.8638 )

 For some more background there was a discussion about the possible
 need for a distinction in NetCDF between cell and corner registration
 here:

 http://mailman.cgd.ucar.edu/pipermail/cf-metadata/2009/002721.html

 In particular this statement is one confirmation I found of the lack
 of cell/corner distinction apart from boundary metadata.

 Coordinate bounds allow edges to be explicit, and you could ask the
 application developer to look
 for those. When not specified, I put the edges halfway between the
 coordinates.

 Cheers, Mike.

 On Tue, Apr 6, 2010 at 3:32 PM, Pinner, Luke
 luke.pin...@environment.gov.au wrote:
  The NetCDF driver computes the pixel centre from lat/lon variables (see
  http://trac.osgeo.org/gdal/ticket/1506).
 
  Is this correct or should it compute pixel upper-left?
 
  Using the example NetCDF dataset attached to the above ticket
  (http://trac.osgeo.org/gdal/attachment/ticket/1506/out.nc), gdalinfo
  reports the min longitude value to be -1.40625 which is not a valid
  longitude (it's 1/2 a pixel less than 0 which is the min value in the
  NetCDF lon variable).
 
  gdalinfo out.nc
  $ gdalinfo -nomd out.nc
  Warning 1: Latitude grid not spaced evenly.
  Seting projection for grid spacing is within 0.1 degrees threshold.
 
  Driver: netCDF/Network Common Data Format
  Files: out.nc
  Size is 128, 64
  Coordinate System is `'
  Origin = (-1.4062500,89.258462312871046)
  Pixel Size = (2.8125000,-2.789326947277220)
  Corner Coordinates:
  Upper Left  (  -1.4062500,  89.2584623)
  Lower Left  (  -1.4062500, -89.2584623)
  Upper Right ( 358.594,  89.258)
  Lower Right ( 358.594, -89.258)
  Center  ( 178.5937500,   0.000)
  Band 1 Block=128x1 Type=Float32, ColorInterp=Undefined
   NoData Value=9.96920996838686905e+36
 
  Regards
  Luke
 
 
  --
  If you have received this transmission in error please notify us
 immediately by return e-mail and delete all copies. If this e-mail or any
 attachments have been sent to you in error, that error does not constitute
 waiver of any confidentiality, privilege or copyright in respect of
 information in the e-mail or attachments.
 
 
 
  Please consider the environment before printing this email.
 
  --
 
  ___
  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] Getting an WMS without XML file.

[Roger André]

Doe this not work for you?

$ python
 from osgeo import gdal
 image = gdal.Open('
http://localhost/cgi-bin/mapserv?map=/var/www/mapfiles/gis_layers/gis_layers.mapREQUEST=GetMapSERVICE=WMSversion=%221.1.122REQUEST=GetMapLAYERS=imagerySTYLES=BBOX=73.617203,18.168884,134.773590,53.554436WIDTH=800HEIGHT=800FORMAT=image/png;%20mode=24bitSRS=epsg:4326
')
 image.RasterXSize
800
 image.RasterYSize
800

Roger
--



On Tue, Apr 6, 2010 at 1:33 AM, Manolo Padron Martinez
manol...@gmail.comwrote:

 Hi:

 I'm trying to use gdal to get images from a WMS server, I got a simple
 code to do it (using GDALOpen and RasterIO methods). But I need to
 make the petition without create previously a XML file, I've been
 looking into documentation but I didn't find a way to do it. Anyone
 knows how to do it?

 Regards from Canary Islands

 Manolo Padron Martinez
 ___
 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] read envisat with gdal

[Roger André]

Can you make a small portion of your data available to us?  If GDAL can read
the data file, then it should be able to read this information as well.

Roger
--

On Fri, Apr 2, 2010 at 5:44 PM, menglong yan yanmengl...@gmail.com wrote:

 Hi, I want to read envisat with gdal, and I can read the data
 successfully. However I can not get the header information such as
 FIRST_MID_LAT(this item could be
 found with  NEST(Next ESA SAR toolbox)) with
 GDALDataset::GetMetadataItem(FIRST_MID_LAT). I have read the souce
 code in /frmts/envisat/, and I am sure it is
 right to call this function. I turn to Anybody who once did this job
 for some advice. Thank you!
 Kind regards.
 ___
 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] read envisat with gdal

[Roger André]

Downloaded a sample file from
http://envisat.esa.int/services/sample_products/aatsr/L2/.  Seems to work
fine in Python.

$ python
 from osgeo import gdal
 image =
gdal.Open('ATS_MET_2PNETB20020729_070738_58362008_00092_02150_0873.N1')
 image.GetMetadata()
{'MPH_STATE_VECTOR_TIME': '29-JUL-2002 08:14:40.282674', snip}
 metadata = image.GetMetadata()
 metadata['SPH_FIRST_MID_LAT']
'+0058766642'
 image.GetMetadataItem('SPH_FIRST_MID_LAT')
'+0058766642'

Can you run gdalinfo on your file and post the output?

Thanks,

Roger
--


On Tue, Apr 27, 2010 at 9:48 AM, Roger André ran...@gmail.com wrote:

 Can you make a small portion of your data available to us?  If GDAL can
 read the data file, then it should be able to read this information as well.

 Roger
 --


 On Fri, Apr 2, 2010 at 5:44 PM, menglong yan yanmengl...@gmail.comwrote:

 Hi, I want to read envisat with gdal, and I can read the data
 successfully. However I can not get the header information such as
 FIRST_MID_LAT(this item could be
 found with  NEST(Next ESA SAR toolbox)) with
 GDALDataset::GetMetadataItem(FIRST_MID_LAT). I have read the souce
 code in /frmts/envisat/, and I am sure it is
 right to call this function. I turn to Anybody who once did this job
 for some advice. Thank you!
 Kind regards.
 ___
 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] Can VRT be used to control display order of images in mosaic?

[Roger André]

Hi All,

I'd like to control the order in which overlapping images are displayed in a
mosaic.  Is it possible to do that by using a VRT, and altering the order in
which the files are defined in it?

The online doc for gdalbuildvrt states, If there is some amount of spatial
overlapping between files, the order may depend on the order they are
inserted in the VRT file, but this behaviour should not be relied on.

Not sure if this disclaimer is for the gdalbuildvrt utility, or describes
the generic behavior of a raster VRT.

Thanks,

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

Re: [gdal-dev] Can VRT be used to control display order of images in mosaic?

[Roger André]

Hi Frank,

Thanks for the reply.  It might be worth keeping the option to control layer
ordering somehow in a mosaic.  The reason I'm exploring this is to see if
it's possible to emulate the ArcGIS 10 closest to center behavior.  I'm
thinking that if I dynamically write the VRT, I can control how images are
overlapped in response to the BBOX extents specified in a WMS request.

Roger
--

On Tue, Jun 22, 2010 at 8:16 AM, Frank Warmerdam warmer...@pobox.comwrote:

 Roger André wrote:

 Hi All,

 I'd like to control the order in which overlapping images are displayed in
 a mosaic.  Is it possible to do that by using a VRT, and altering the order
 in which the files are defined in it?

 The online doc for gdalbuildvrt states, If there is some amount of
 spatial overlapping between files, the order may depend on the order they
 are inserted in the VRT file, but this behaviour should not be relied on.

 Not sure if this disclaimer is for the gdalbuildvrt utility, or describes
 the generic behavior of a raster VRT.


 Roger,

 I believe the statement applies to the VRT format itself.  I have looked
 at the code, and it always applies them linearly based on the order they
 appear in the .VRT list of sources.  I think you could depend on it for
 now
 though in a future version we might do something more fancy that would
 result in alternate ordering.

 Best regards,
 --

 ---+--
 I set the clouds in motion - turn up   | Frank Warmerdam,
 warmer...@pobox.com
 light and sound - activate the windows | 
 http://pobox.com/~warmerdamhttp://pobox.com/%7Ewarmerdam
 and watch the world go round - Rush| Geospatial Programmer for Rent


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

[gdal-dev] Changing precision of coordinates in Python OGR API

[Roger André]

Hi All,

I have an assignment that requires me to truncate the precision of all my
polygon coordinates to 6 decimal places.  Is there an easy way to do this
while I have the feature opened with OGR and can access its ogr.Geometry
directly?  If not, can someone recommend a good way to iterate through each
coordinate in a feature in such a way that I can treat it as a string and
format it?

Thanks,

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

Re: [gdal-dev] Changing precision of coordinates in Python OGR API

[Roger André]

Thanks Chaitanya,

I like the look of that second option.  If I understand correctly what Even
says, I should then be able to call the *ExportToWKT() *method with a
precision parameter?

Roger
--

On Tue, Aug 16, 2011 at 10:36 PM, Chaitanya kumar CH chaitanya...@gmail.com
 wrote:

 Roger,

 OGR doesn't support control of coordinate precision. You need to do this
 yourselves.
 You can get the WKT of the geometry and and edit the string. If you are
 willing to recompile GDAL/OGR, you can modify the code that creates the WKT.
 Even wrote how to do this in this mailing list [1].

 [1]: http://lists.osgeo.org/pipermail/gdal-dev/2011-April/028306.html

 On Wed, Aug 17, 2011 at 5:21 AM, Roger André ran...@gmail.com wrote:

 Hi All,

 I have an assignment that requires me to truncate the precision of all my
 polygon coordinates to 6 decimal places.  Is there an easy way to do this
 while I have the feature opened with OGR and can access its ogr.Geometry
 directly?  If not, can someone recommend a good way to iterate through each
 coordinate in a feature in such a way that I can treat it as a string and
 format it?

 Thanks,

 Roger

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




 --
 Best regards,
 Chaitanya kumar CH.

 +91-9494447584
 17.2416N 80.1426E

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

Re: [gdal-dev] Changing precision of coordinates in Python OGR API

[Roger André]

Works brilliantly, thanks again.

- Set env var:
export OGR_WKT_PRECISION=4

- In Python interpreter:
from osgeo import ogr
datasource = ogr.Open('World_Admin_2.shp')
layer = datasource.GetLayer(0)
feat = layer.GetNextFeature()
geom = feat.geometry()
wkt = geom.ExportToWkt()

was:
151.85768502434 -32.613605699642989,151.857672992 -32.613692002142983))
now:
151.8577 -32.6136,151.8577 -32.6137))

Roger
---

On Tue, Aug 16, 2011 at 10:53 PM, Roger André ran...@gmail.com wrote:

 Thanks Chaitanya,

 I like the look of that second option.  If I understand correctly what Even
 says, I should then be able to call the *ExportToWKT() *method with a
 precision parameter?

 Roger
 --

 On Tue, Aug 16, 2011 at 10:36 PM, Chaitanya kumar CH 
 chaitanya...@gmail.com wrote:

 Roger,

 OGR doesn't support control of coordinate precision. You need to do this
 yourselves.
 You can get the WKT of the geometry and and edit the string. If you are
 willing to recompile GDAL/OGR, you can modify the code that creates the WKT.
 Even wrote how to do this in this mailing list [1].

 [1]: http://lists.osgeo.org/pipermail/gdal-dev/2011-April/028306.html

 On Wed, Aug 17, 2011 at 5:21 AM, Roger André ran...@gmail.com wrote:

 Hi All,

 I have an assignment that requires me to truncate the precision of all my
 polygon coordinates to 6 decimal places.  Is there an easy way to do this
 while I have the feature opened with OGR and can access its ogr.Geometry
 directly?  If not, can someone recommend a good way to iterate through each
 coordinate in a feature in such a way that I can treat it as a string and
 format it?

 Thanks,

 Roger

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




 --
 Best regards,
 Chaitanya kumar CH.

 +91-9494447584
 17.2416N 80.1426E



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

Re: [gdal-dev] KMLSuperOverlay in Google Earth 6.1.0.5001

[Roger André]

What is the output of gdalinfo on one of the PNG's in the superoverlay?



On Thu, Nov 17, 2011 at 1:36 PM, Roland Duhaime roland.duha...@gmail.comwrote:


 Has anyone else experienced the latest version of Google Earth
 (6.1.0.5001) or the G.E. API not loading PNG formatted KMLSUPEROVERLAYs
 created using gdal_translate?  The same KMZ files worked in the past.

 Thanks,
 Roland


 ___
 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] What's up with the 8859 encoding?

[Roger André]

Even,

Thank you for the quick reply.   And it appears that using a .CPG file does
not make a difference in this.  Is that correct?

Roger

On Fri, Nov 7, 2014 at 3:00 PM, Even Rouault even.roua...@spatialys.com
wrote:

 Le vendredi 07 novembre 2014 23:53:09, Roger André a écrit :
  Hello List,
 
  I am very curious why after installing GDAL 1.11, I now need to use the
  -lco ENCODING=UTF-8 option when using ogr2ogr on any of my already UTF-8
  encoded shapefiles?  Failure to do so results in this error,
 
  Warning 1: One or several characters couldn't be converted correctly from
  UTF-8 to ISO-8859-1.

 Roger,

 The default encoding of shapefiles is ISO-8859-1. So if you don't specify
 -lco
 ENCODING=UTF-8, the shapefile driver will attempt to recode from UTF-8 to
 ISO-8859-1. This behaviour has been introduced in GDAL 1.9. Previously no
 transcoding was done, and there could be mismatch between the actual
 encoding
 of the content and the declared encoding of the DBF.

 Even

 --
 Spatialys - Geospatial professional services
 http://www.spatialys.com

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

Re: [gdal-dev] What's up with the 8859 encoding?

[Roger André]

I see.  So the output will default to ISO-8859-1 - regardless of the input
.CPG - if the ENCODING flag isn't used.

On Fri, Nov 7, 2014 at 3:08 PM, Even Rouault even.roua...@spatialys.com
wrote:

 Le samedi 08 novembre 2014 00:03:22, Roger André a écrit :
  Even,
 
  Thank you for the quick reply.   And it appears that using a .CPG file
 does
  not make a difference in this.  Is that correct?

 You have to realize that the input and output side of ogr2ogr are
 completely
 isolated from each other.

 So if your input shapefile has a .cpg,  the shapefile driver will use it to
 recode from the input dbf encoding mentionned in the .cpf to UTF-8 which is
 the pivot encoding of OGR (if the .cpg indicates UTF-8 then no recoding
 occurs)
 If the output is a shapefile and you specify -lco ENCODING=UTF-8, as the
 pivot
 encoding of OGR is UTF-8, no recoding will occur and an output .cpg file
 indicating UTF-8 will be generated.

 
  Roger
 
  On Fri, Nov 7, 2014 at 3:00 PM, Even Rouault even.roua...@spatialys.com
 
 
  wrote:
   Le vendredi 07 novembre 2014 23:53:09, Roger André a écrit :
Hello List,
   
I am very curious why after installing GDAL 1.11, I now need to use
 the
-lco ENCODING=UTF-8 option when using ogr2ogr on any of my already
UTF-8 encoded shapefiles?  Failure to do so results in this error,
   
Warning 1: One or several characters couldn't be converted correctly
from UTF-8 to ISO-8859-1.
  
   Roger,
  
   The default encoding of shapefiles is ISO-8859-1. So if you don't
 specify
   -lco
   ENCODING=UTF-8, the shapefile driver will attempt to recode from UTF-8
 to
   ISO-8859-1. This behaviour has been introduced in GDAL 1.9. Previously
 no
   transcoding was done, and there could be mismatch between the actual
   encoding
   of the content and the declared encoding of the DBF.
  
   Even
  
   --
   Spatialys - Geospatial professional services
   http://www.spatialys.com

 --
 Spatialys - Geospatial professional services
 http://www.spatialys.com

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

[gdal-dev] What's up with the 8859 encoding?

[Roger André]

Hello List,

I am very curious why after installing GDAL 1.11, I now need to use the
-lco ENCODING=UTF-8 option when using ogr2ogr on any of my already UTF-8
encoded shapefiles?  Failure to do so results in this error,

Warning 1: One or several characters couldn't be converted correctly from
UTF-8 to ISO-8859-1.

What gives?

Thanks,

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

Re: [gdal-dev] Aligning Raster Over the Polygon

[Roger André]

Have you tried using a World File?  http://en.wikipedia.org/wiki/World_file

Roger

On Tue, Nov 18, 2014 at 10:29 AM, Simen Langseth simlan...@gmail.com
wrote:

 Dear GDAL Users:

 I have a raster image which have to be aligned over the polygon.
 Both the raster and polygon have the same coordinate system.

 The figure shows the raster image which has to be fitted over the RED
 polygon.


 https://drive.google.com/file/d/0B2RqG9tSAAIUc3J5ZDN6V01IcG8/view?usp=sharing

 How can I do it guys?

 Thanks.

 Simen

 ___
 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] georeferencing raster data using vector polygon

[Roger André]

Hi Simen,

Ok, here is how to implement the suggestion I gave, which was to create a
new worldfile and shift the image so that it fit the footprint vector
better.

1. Use gdal_edit.py with the -unsetgt option to remove the georeferencing
from the .tif
2. Then copy the attached .tfw file into the same dir as the .tif and view
it in relation to the KML footprint.
3.  You can tweak the alignment further by modifiying the bottom 2
coordinate values in the .tfw file.

Why do I advocate this method over warping the image using GCP's or QGis?
Because those are lossy processes.  In other words, the warping will
force a re-write of the pixel values using some sort of interpolation.  You
might still want to do that at the very end, like if you need  perfect
point-to-point matching against the footprint corners, but I still would
shift the image first.

Hope this helps,

Roger

On Wed, Nov 19, 2014 at 1:16 AM, Simen Langseth simlan...@gmail.com wrote:

  I wanted to georeference the raster image in UTM coordinates into the
 boundry of red polygon with geographic coordinate system.


 https://drive.google.com/file/d/0B2RqG9tSAAIUc3J5ZDN6V01IcG8/view?usp=sharing

 I tried as follows:

 gdalwarp,-t_srs wgs84 -srcnodata 0 -dstnodata 0 -r near -te -96.75234
 54.93479 -92.90071 56.91041 i

 What parameters should I provide for that?

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



p031r021.tfw
Description: Binary data
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Setting encoding using Python API

[Roger André]

Hi All,

Curious how I go about specifying the equivalent of -lco ENCODING=UTF-8
when using the Python API to write out shapefiles with OGR?

Thanks,

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

Re: [gdal-dev] Setting encoding using Python API

[Roger André]

Awesome!  I hate to ask [another] dumb question, but how can I peruse that
autotest suite?  Under some sort of test dir in the install?

Roger

On Fri, Jan 30, 2015 at 1:39 AM, Even Rouault even.roua...@spatialys.com
wrote:

 Le vendredi 30 janvier 2015 10:37:30, Roger André a écrit :
  Hi All,
 
  Curious how I go about specifying the equivalent of -lco ENCODING=UTF-8
  when using the Python API to write out shapefiles with OGR?

 Roger,

 sure: taken from autotest suite:

 shape_ds.CreateLayer('test_utf_cut', geom_type = ogr.wkbPoint, \
 options = ['ENCODING=UTF-8'])


 Even

 --
 Spatialys - Geospatial professional services
 http://www.spatialys.com

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

Re: [gdal-dev] [Proj] Pixel Shifts When Warping

[Roger André]

The GDAL data model does specifically call out the top-left corner of the
top-left pixel as the origin. Not sure if that's helpful, but I remember
having to look this up before.

http://www.gdal.org/gdal_datamodel.html


On Mon, May 18, 2015 at 8:31 AM, Yuta Sato yutaxs...@gmail.com wrote:

 Dear all:

 Thanks for sharing your ideas and experiences.

 As I mentioned in my question, it seems that MRT result is shifting right
 while the GDAL warp result is correct. Do you agree?

 The exact upper left corner of the original HDF-EOS file is exactly as the
 GDAL warpped result.

 The HDF file uses top left corner of top left pixel, while the MRT seems
 to use central of top left pixel (MRT userguide page 10 under corner
 coordinates section), am I right?


 https://lpdaac.usgs.gov/sites/default/files/public/mrt41_usermanual_032811.pdf

 If I am right, fortunately the GDAL would win!





 On Mon, May 18, 2015 at 11:24 PM, Markus Neteler nete...@osgeo.org
 wrote:

 On Sun, May 17, 2015 at 1:32 PM, Even Rouault
 even.roua...@spatialys.com wrote:
  Le dimanche 17 mai 2015 12:54:44, Markus Neteler a écrit :
  On Sun, May 17, 2015 at 6:44 AM, Yuta Sato yutaxs...@gmail.com
 wrote:
   Dear GDAL and Proj Developers and Users:
  
   I am really frustrated when the MODIS HDF file converted into GeoTiff
   with GCS projection never match with the same MODIS file  converted
 from
   MODIS Reprojection tool.
 
  We had this issue in the past, too, for our EuroLST project
  (gap-filling of all LST maps for Europe, see
 http://gis.cri.fmach.it/eurolst/).
  Solved since then, see below.
  ...
 
   Do you have any idea on which is more precise: gdalwarp or MRT? Or
 there
   is any way in gdalwarp to get the same image as MRT gives?

 An addition:

 For our EuroLST project we use EU LAEA as target projection. In our
 tests (we made a lot) the MRT results came out shifted by 0.5 -1.0
 pixel since MRT does not support a datum or ellipsoid for LAEA (LA in
 MRT), instead it supports only a sphere when reprojecting from MODIS
 sinusoidal.

 Using recent GDAL, gdalwarp does the reprojection from MODIS
 sinusoidal to EU LAEA correctly.

  You may take a look at
  http://pymodis.fem-environment.eu/
 
  pyModis = 1.0 comes with fast GDAL support (way faster than MRT).
 
  The details of the gdalwarp usage you find in the source code here:
 
 https://github.com/lucadelu/pyModis/blob/master/pymodis/convertmodis_gdal.py
 
  Markus,
 
  for the sake of my curiosity, could you explain what specific stuff it
 does
  compared to gdalwarp ? I quickly skimmed through the file and only saw
 regular
  use of the GDAL warping API. But perhaps I missed something.

 pyMODIS is a set of scripts and a library to perform the download of
 large numbers of MODIS HDF/XML files (think cronjob), it creates a
 mosaic of several tiles, it extracts specific information from
 bit-encoded MODIS quality assessment layers of different product types
 etc.

 Essentially a ready to use product which happily uses GDAL as a backend.

 Markus

 --
 http://courses.neteler.org
 http://gis.cri.fmach.it/neteler/
 http://grass.osgeo.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

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

[gdal-dev] Problem appending non-geometry tables in PostGis

[Roger André]

Hi All,

I'm having some trouble using ogr2ogr to do batch uploads to a Postgis DB.
It appears that on tables which don't contain geometry, CreateLayer is
failing and not allowing data to be appended to an existing table.  Here is
the command that is being used to load each feature class:

ogr2ogr -f PostgreSQL -lco ENCODING=UTF-8 -append -update PG:host=foo
dbname=bar user=crak password=head ./nokia_here_asia/file.gdb CityPOINames

Below is my logging output.  The first set of data shows the asia data set
being loaded, which creates the PostGIS tables. The next set shows what
happens with the europe data set.  Although warnings are given for
MapAdminArea, MapAdminLink and CityPOI layers, they are sucessfully
appended.  However, CityPOINames, MapAreaNames, and MapNameTrans all fail
with the CreateLayer error.  Those tables are different in that they lack
geometry.

Extracting from ./nokia_here_asia/file.gdb
Loading MapAdminArea feature_class
Loading MapAdminLink feature_class
Loading CityPOI feature_class
Loading CityPOINames feature_class
Loading MapAreaNames feature_class
Loading MapNameTrans feature_class

Extracting from ./nokia_here_europe/file.gdb
Loading MapAdminArea feature_class
WARNING: Layer creation options ignored since an existing layer is
 being appended to.
Loading MapAdminLink feature_class
WARNING: Layer creation options ignored since an existing layer is
 being appended to.
Loading CityPOI feature_class
WARNING: Layer creation options ignored since an existing layer is
 being appended to.
Loading CityPOINames feature_class
ERROR 1: Layer citypoinames already exists, CreateLayer failed.
Use the layer creation option OVERWRITE=YES to replace it.
ERROR 1: Terminating translation prematurely after failed
translation of layer CityPOINames (use -skipfailures to skip errors)

Loading MapAreaNames feature_class
ERROR 1: Layer mapareanames already exists, CreateLayer failed.
Use the layer creation option OVERWRITE=YES to replace it.
ERROR 1: Terminating translation prematurely after failed
translation of layer MapAreaNames (use -skipfailures to skip errors)

Loading MapNameTrans feature_class
ERROR 1: Layer mapnametrans already exists, CreateLayer failed.
Use the layer creation option OVERWRITE=YES to replace it.
ERROR 1: Terminating translation prematurely after failed
translation of layer MapNameTrans (use -skipfailures to skip errors)


I could create some giant CSV files and load these in one shot, but I'd
prefer it if I could handle the entire load using ogr2ogr.  Has anyone else
encountered this, or have suggestions for dealing with the problem?

Thanks,

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

Re: [gdal-dev] Problem appending non-geometry tables in PostGis

[Roger André]

Even,

That works, thanks.

Roger

On Fri, Jun 5, 2015 at 12:52 AM, Even Rouault even.roua...@spatialys.com
wrote:

 Le vendredi 05 juin 2015 02:19:06, Roger André a écrit :
  Hi All,
 
  I'm having some trouble using ogr2ogr to do batch uploads to a Postgis
 DB.
  It appears that on tables which don't contain geometry, CreateLayer is
  failing and not allowing data to be appended to an existing table.  Here
 is
  the command that is being used to load each feature class:
 
  ogr2ogr -f PostgreSQL -lco ENCODING=UTF-8 -append -update PG:host=foo
  dbname=bar user=crak password=head ./nokia_here_asia/file.gdb
 CityPOINames

 Roger,

 Add --config PG_LIST_ALL_TABLES YES

 Even

 --
 Spatialys - Geospatial professional services
 http://www.spatialys.com

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

Re: [gdal-dev] Merge all features in polygon shapefile with gdal

[Roger André]

Hey Steven,

Are you trying to do this via a spatially enabled database?  In PostGIS,
the dissolve syntax looks like this:

SELECT dissolve_field, ST_UNION(geometry) FROM table GROUP BY
dissolve_field;

Roger


On Mon, Jul 6, 2015 at 4:13 PM, Steven Lutz ste...@ruggedrocksoffroad.com
wrote:

  I've found some examples of how to dissolve features in a polygon using
 gdal but I'm not familiar with how to implement them.
 here is an example of what I've found, but I'm having trouble modifying it

 this is what I've come up with:

 ogr2ogr output.shp
   input.
 shp -dialect sqlite -sql
 SELECT ST_Union(geometry), dissolve_field FROM input GROUP BY dissolve_field


 it looks to me that it is trying to do this operation based on specific
 attribute entries but I'm looking to ignore the attribute entrie and make
 all features into one.
 This is what I've tried:
 the 32768.shp has an attribute named BLOCK and it's still not working..
 any ideas? ogr2ogr -f ESRI Shapefile output.shp 32768.shp -dialect
 sqlite -sql select ST_union(block) from '32768' Warning 6:
 Normalized/laundered field name: 'ST_union(block)' to 'ST_union(b

 I am confused as to what I should be putting in the parenthasis of
 ST_Union()

 Please advise.
 Thank you


 --
  Steven Lutz
 President
 Rugged Rocks, LLC
 Sales  Support: (909)547-4651
 Nissan 4x4 Parts http://ruggedrocks.com / Off Road GPS Units  Maps
 http://ruggedroutes.com

 ___
 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] Weird OpenFileGDB error

[Roger André]

Hi All,

I've been trying unsuccessfully to extract a feature class from a large
file geodatabase.  While many of the other feature classes in it are read
properly, on this one I get the following error:

$ ogrinfo -summary file.gdb Streets
Had to open data source read-only.
INFO: Open of `file.gdb'
  using driver `OpenFileGDB' successful.

Layer name: Streets
Geometry: Multi Line String
ERROR 1: Error occured in filegdbtable.cpp at line 768

Line 768 is part of this block of code:

 if( byFieldType > FGFT_XML)
{
CPLDebug("OpenFileGDB", "Unhandled field type : %d",
byFieldType);
returnError();   < line 768
}

Can any of you provide me some suggestions for how I can trace the specific
cause of this problem more precisely?  Unfortunately, this is commercial
data which I can't share with you, but I'd be happy to try to diagnose the
problem more thoroughly.

Thanks,

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

Re: [gdal-dev] How to write WKB as string?

[Roger André]

Ahh, thanks!  It helps to frame the question correctly.  I see now that
what I am trying to do is create the hex-encoded form of WKB.

On Wed, Oct 14, 2015 at 10:24 AM, Even Rouault <even.roua...@spatialys.com>
wrote:

> Le mercredi 14 octobre 2015 19:18:20, Roger André a écrit :
> >   Yes Even, exportToWkb() was something I tried, but I thought
> there
> > might be an easy way to get the ASCII form of the WKB, rather than the
> > straight binary that I got when I wrote out to a text file.
>
> Ah I missed that part.
>
> The following should help:
>
> http://www.gdal.org/cpl__string_8h.html:
>
> char *  CPLBinaryToHex (int nBytes, const GByte *pabyData)
> CPL_WARN_UNUSED_RESULT
> Binary to hexadecimal translation. More...
>
> GByte * CPLHexToBinary (const char *pszHex, int *pnBytes)
> CPL_WARN_UNUSED_RESULT
>
> Even more compact with base64 encoding:
>
> char *  CPLBase64Encode (int nBytes, const GByte *pabyData)
> CPL_WARN_UNUSED_RESULT
>
> int CPLBase64DecodeInPlace (GByte *pszBase64)
>
> >
> > Norman, thanks for the links. I'll go peruse them.
> >
> > On Wed, Oct 14, 2015 at 10:10 AM, Norman Vine <n...@cape.com> wrote:
> > > On Oct 14, 2015, at 12:59 PM, Roger André <ran...@gmail.com> wrote:
> > >
> > > Hi All,
> > >
> > > I have a tool which currently exports geometry using the ExportToWkt()
> > > function.  It works great, but I'd like to try reducing the size of my
> > > output files by using a format like the string representation of WKB
> that
> > > shp2pgsql produsces when creating SQL insert files.  Is there a export
> > > function available in GDAL that I could use for this?
> > >
> > > I'm looking at http://www.gdal.org/classOGRGeometry.html currently,
> but
> > > have not yet spotted what I'm looking for.
> > >
> > > Thanks,
> > >
> > >
> > > see
> > > http://www.gdal.org/ogr__api_8h.html#ad5c94ae76c09774dba8725c24daeefd6
> > > or
> > >
> http://www.gdal.org/classOGRGeometry.html#a9a3ad2f9b55285be500ca6ac5dc75d
> > > 1d
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] How to write WKB as string?

[Roger André]

Hi All,

I have a tool which currently exports geometry using the ExportToWkt()
function.  It works great, but I'd like to try reducing the size of my
output files by using a format like the string representation of WKB that
shp2pgsql produsces when creating SQL insert files.  Is there a export
function available in GDAL that I could use for this?

I'm looking at http://www.gdal.org/classOGRGeometry.html currently, but
have not yet spotted what I'm looking for.

Thanks,

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

Re: [gdal-dev] How to write WKB as string?

[Roger André]

  Yes Even, exportToWkb() was something I tried, but I thought there
might be an easy way to get the ASCII form of the WKB, rather than the
straight binary that I got when I wrote out to a text file.

Norman, thanks for the links. I'll go peruse them.

On Wed, Oct 14, 2015 at 10:10 AM, Norman Vine <n...@cape.com> wrote:

>
> On Oct 14, 2015, at 12:59 PM, Roger André <ran...@gmail.com> wrote:
>
> Hi All,
>
> I have a tool which currently exports geometry using the ExportToWkt()
> function.  It works great, but I'd like to try reducing the size of my
> output files by using a format like the string representation of WKB that
> shp2pgsql produsces when creating SQL insert files.  Is there a export
> function available in GDAL that I could use for this?
>
> I'm looking at http://www.gdal.org/classOGRGeometry.html currently, but
> have not yet spotted what I'm looking for.
>
> Thanks,
>
>
> see
> http://www.gdal.org/ogr__api_8h.html#ad5c94ae76c09774dba8725c24daeefd6
> or
> http://www.gdal.org/classOGRGeometry.html#a9a3ad2f9b55285be500ca6ac5dc75d1d
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Python feature.SetField() fails with unicode strings

[Roger André]

Good morning,

I've been working on a simple script to load GeoJSON into PostGIS.  I'm
encountering an error that has me scratching my head, it seems that unicode
strings cannot be used with the feature.SetField() function in Python.  Is
that really true?

I've tested this by re-writing the variable inline from my test record and
can confirm that when it is type  it works, and when it is type
 it fails with the following error.

 *NotImplementedError: Wrong number of arguments for overloaded function
'Feature_SetField'.*

I foind an old ticket that talks about this,
https://trac.osgeo.org/gdal/ticket/4608, but I'm not having any luck
getting this to work.

My GDAL version info is GDAL 1.11.0, released 2014/04/16.

To reiterate, I'm simply parsing a field from a UTF-8 encoded JSON record
and trying to push it into PostGIS.

Thank you for your time.

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

Re: [gdal-dev] Python feature.SetField() fails with unicode strings

[Roger André]

Hi Even,

Thank you for this information.  I'll see about getting the latest version
of GDAL installed on my system.  I tried using the string.encode('utf-8')
workaround and continue to get the same error.

Roger

On Mon, Jul 24, 2017 at 8:46 AM, Even Rouault <even.roua...@spatialys.com>
wrote:

> On lundi 24 juillet 2017 08:43:02 CEST Roger André wrote:
>
> > Good morning,
>
> >
>
> > I've been working on a simple script to load GeoJSON into PostGIS. I'm
>
> > encountering an error that has me scratching my head, it seems that
> unicode
>
> > strings cannot be used with the feature.SetField() function in Python. Is
>
> > that really true?
>
> >
>
> > I've tested this by re-writing the variable inline from my test record
> and
>
> > can confirm that when it is type  it works, and when it is type
>
> >  it fails with the following error.
>
> >
>
> > *NotImplementedError: Wrong number of arguments for overloaded function
>
> > 'Feature_SetField'.*
>
> >
>
> > I foind an old ticket that talks about this,
>
> > https://trac.osgeo.org/gdal/ticket/4608, but I'm not having any luck
>
> > getting this to work.
>
>
>
> This ticket was fixed in GDAL 2.0. There might have been other fixes
> since. So basically you should try upgrading. Or as a workaround use
> your_string.encode('utf-8') to transform from Unicode to UTF-8 encoded
> string.
>
>
>
> Even
>
>
>
> --
>
> Spatialys - Geospatial professional services
>
> http://www.spatialys.com
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Large number of "Unable to translate coordinate system" in epsg file

[Roger André]

OK, thank you.

On Tue, Apr 17, 2018, 4:53 AM Kristian Evers <kr...@sdfe.dk> wrote:

> I've added a ticket for the addition for the Michigan LCC:
> https://github.com/OSGeo/proj.4/issues/940
>
> It is a rather trivial change so a good issue for a first time PROJ
> contributor.
>
> /Kristian
>
> -Oprindelig meddelelse-
> Fra: gdal-dev [mailto:gdal-dev-boun...@lists.osgeo.org] På vegne af Even
> Rouault
> Sendt: 17. april 2018 10:49
> Til: gdal-dev@lists.osgeo.org
> Emne: Re: [gdal-dev] Large number of "Unable to translate coordinate
> system" in epsg file
>
> On lundi 16 avril 2018 10:46:18 CEST Roger André wrote:
> > Hi All,
> >
> > Slightly off-topic question, but I'm curious why I'm seeing so many
> entries
> > that look like this in my epsg file:
> >
> > # NAD27 / Michigan North
> > # Unable to translate coordinate system EPSG:6200 into PROJ.4 format.
> >
> > Is this normal, or the sign of a flawed install of dependencies, or
> > something else?
>
> Yes this is "normal".
>
> EPSG:6200 / "NAD27 / Michigan North" uses a LCC_2SP(Michigan) projection
> (EPSG:1051) which is not supported by GDAL / proj :
> http://lists.maptools.org/pipermail/proj/2015-April/007108.html
>
> Even
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

[gdal-dev] Large number of "Unable to translate coordinate system" in epsg file

[Roger André]

Hi All,

Slightly off-topic question, but I'm curious why I'm seeing so many entries
that look like this in my epsg file:

# NAD27 / Michigan North
# Unable to translate coordinate system EPSG:6200 into PROJ.4 format.

Is this normal, or the sign of a flawed install of dependencies, or
something else?

Thanks,

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

Re: [gdal-dev] Import Shapefile into SQL Server on Ubuntu

[Roger André]

See if this helps -
https://docs.microsoft.com/en-us/sql/connect/odbc/linux-mac/installing-the-microsoft-odbc-driver-for-sql-server?view=sql-server-2017
.

On Mon, Oct 22, 2018 at 8:25 AM Peter Marlow 
wrote:

> Hi,
>
>
>
> I’m trying to import a Shapefile into an SQL Server database using ogr2ogr
> on Ubuntu. The command I’m running looks like:
>
>
>
> *ogr2ogr -f MSSQLSpatial
> "MSSQL:server=localhost;database=[database-name];username=sa;password=*;trusted_connection=yes"
> [path-to-shape-file]*
>
>
>
> It returns an error stating:
>
>
>
> *ERROR 1: Unable to find driver `MSSQLSpatial'.*
>
>
>
> Is it possible to get the driver on linux? Or use an alternative driver?
> I’ve had a google and can’t seem to see any workarounds.
>
>
>
> The docs here (https://www.gdal.org/drv_mssqlspatial.html) talk about
> specifying a *Driver* parameter that could be a custom SQL Server driver,
> however this is in the connection string and not as a parameter for the
> ogr2ogr command.
>
>
>
> Thanks,
> Pete
>
>
>
> 
>
>
> SCISYS UK Limited. Registered in England and Wales No. 4373530.
> Registered Office: Methuen Park, Chippenham, Wiltshire SN14 0GB, UK.
>
> Before printing, please think about the environment.
>
> 
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev