[gdal-dev] profiles in gdal2tiles.py
Hi All, I'm interested in experimenting with creating tile sets in various coordinate systems (other than the usual EPSG:3857 and EPSG:4326). The good news is that it looks like as of GDAL v3.2 gdal2tiles.py supports a "profile" option for defining the tile set including the coordinate system. The profile is in JSON with all the detail gdal2tiles would require to create tiles (example: https://github.com/OSGeo/gdal/blob/master/gdal/data/tms_MapML_APSTILE.json) I suppose I could create a profile JSON file by hand (or do some coding) but I thought I would ask if there was an existing tool? Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Problem in ogr converting polygon coverages
You're awesome Even! Thanks! Brent From: "Even Rouault" Sent: Thursday, June 29, 2017 2:29 PM To: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] Problem in ogr converting polygon coverages On jeudi 29 juin 2017 21:28:19 CEST Even Rouault wrote: > Brent, > > > I'm having trouble with OGR converting ESRI polygon coverages to > > shapefiles. Basically the newer versions of ogr ignore the attributes > > on the polygons and show only the ArcIds attribute. Looks like a bug > > possibly introduced during RFC 50 > > (https://trac.osgeo.org/gdal/wiki/rfc50_ogr_field_subtype). > > I don't think this is related. More likely a regression caused by a code > cleanup/hardening. > > If there is a option to the command line to fix this, please let me > > > > know, otherwise I'll file a bug report. > > Please file a bug report Actually I've just done it with https://trac.osgeo.org/gdal/ticket/6950 and committed the fix. 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] Problem in ogr converting polygon coverages
Hi all, I'm having trouble with OGR converting ESRI polygon coverages to shapefiles. Basically the newer versions of ogr ignore the attributes on the polygons and show only the ArcIds attribute. Looks like a bug possibly introduced during RFC 50 (https://trac.osgeo.org/gdal/wiki/rfc50_ogr_field_subtype). If there is a option to the command line to fix this, please let me know, otherwise I'll file a bug report. Details: Ogrinfo in GDAL version 1.11.0 works fine (note all the attributes after ArcIds): C:\Projects>C:\Programs\gdal-1-11-0-mapserver-6-4-1\bin\ogrinfo.exe --version GDAL 1.11.0, released 2014/04/16 C:\Projects>C:\Programs\gdal-1-11-0-mapserver-6-4-1\bin\ogrinfo -geom=no ea_20090105 PAL | more Had to open data source read-only. INFO: Open of `ea_20090105' using driver `AVCBin' successful. Layer name: PAL Geometry: Polygon Feature Count: 467 Extent: (-195031.966405, 2723260.836220) - (1970740.118994, 5186575.947505) Layer SRS WKT: : ArcIds: IntegerList (0.0) AREA: Real (18.5) PERIMETER: Real (18.5) EA_20090105#: Integer (5.0) EA_20090105-ID: Integer (5.0) A_LEGEND: String (11.0) REGION: String (2.0) DATE_CARTE: String (10.0) SOURCE: String (5.0) MOD: String (1.0) EGG_ID: Integer (5.0) PNT_TYPE: Integer (3.0) EGG_NAME: String (2.0) EGG_SCALE: Integer (5.0) EGG_ATTR: String (49.0) USER_ATTR: String (49.0) ROTATION: Integer (3.0) E_CT: String (2.0) E_CA: String (1.0) E_CB: String (1.0) : while GDAL version 2.1.3 just shows the ArcIds attribute: C:\Projects>ogrinfo.exe --version GDAL 2.1.3, released 2017/20/01 C:\Projects>ogrinfo -geom=no ea_20090105 PAL | more Had to open data source read-only. INFO: Open of `ea_20090105' using driver `AVCBin' successful. Layer name: PAL Geometry: Polygon Feature Count: 467 Extent: (-195031.966405, 2843291.569215) - (1870695.841972, 4961994.591590) Layer SRS WKT: : ArcIds: IntegerList (0.0) OGRFeature(PAL):2 ArcIds (IntegerList) = (146:883,-893,-889,-868,-890,-898,-886,-885,-880,-875,-854,-856,-841,-840,...) OGRFeature(PAL):3 ArcIds (IntegerList) = (5:1,12,14,15,2) == To try it out, download the e00 file: http://ice-glaces.ec.gc.ca//www_archive/AOI_11/Coverages/rgc_a11_20090105_CEXPREA.e00 and use "avcimport" to convert it to a coverage avcimport gc_a11_20090105_CEXPREA.e00 ea_20090105 then use the ogrinfo command shown above. -- Best Regards, Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] reading MS Access .accdb
Hi all, I've been trying to get ogrinfo to read some MS Access database (.accdb) files so I can link them to some shapefiles. After some digging, I've found that to get ogrinfo to read them, I need to: 1. Rename the .accdb to .mdb 2. Install the MS Access Database Engine Redistributable, such as https://www.microsoft.com/en-us/download/details.aspx?id=13255 3. Specify the driver string for the 64-bit driver (I'm using Windows 7 64-bit and OSGeo4W 64-bit): --config MDB_DRIVER_TEMPLATE "DRIVER=Microsoft Access Driver (*.mdb, *.accdb);DBQ=%s" 4. Use an explicit path to the file (note the .\ or the full path is required): ogrinfo.exe --config MDB_DRIVER_TEMPLATE "DRIVER=Microsoft Access Driver (*.mdb, *.accdb);DBQ=%s" .\t_408088.mdb just setting the current dir and specifying "t_408088.mdb" does not work. I consider 1. and 4. above to be functionality issues/bugs (and 3. is a documentation opportunity). I see that there was a related issue (and a fix) to use the 64-bit driver https://trac.osgeo.org/gdal/ticket/5594 Any comments before I file my bug reports? -- Best Regards, Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] PDF and EXTRA_STREAM
I've been experimenting with generating Geospatial PDFs with GDAL and would like to wrtie some extra text and lines to the PDF file. The PDF drive has a Create Option of EXTRA_STREAM which looks pretty good. I can add text, lines and rectangles by doing: -co EXTRA_STREAM="BT 0 0 255 rg /FTimesRoman 50 Tf 1 0 0 1 1 1 Tm 100 200 Td (Hello World!) Tj ET 80 260 m 80 50 l h S 0.9 0.5 0.0 rg 400 50 300 200 re f" But I think the length of the stream could get very long. Any interest in adding another Create Option like EXTRA_STREAM_FILE=my_stuff.txt so I could put all the extra in a separate file? -- Best Regards, Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Feature request: gauss and other interpolations in gdalwarp
Jan, I can't speak the development request, but you may be able to do some sharpening using the object in a VRT. I've used -0.111 -0.111 -0.111 -0.111 2 -0.111 -0.111 -0.111 -0.111 to apply sharpening to satellite imagery. And a few years ago I found that using wavelet compression on scanned topographic maps when downsampling made the text more readable (preserves edges more). I used ECW back then, but I suspect JPEG2000 may do it as well... Best Regards, Brent Fraser On 3/28/2013 5:27 AM, Jan Hartmann wrote: Perhaps I should clarify a bit what I meant, I haven't had any reactions until now, positive or negative, and it is important for me. I use Gdalwarp and gdaladdo extensively for goereferencing and tiling large historical maps serieses (raster scans). To display them efficiently, I need to create layers at different scale levels, e.g. if the original maps have a pixel size of 10 meters, I need to resample them to rasters with pixel sizes of 20, 40, 80 , 160, 320 and 640 meters, and tile all those maps appropriately. So the original maps have to be resampled quite drastically. In Gdalwarp there is no adequeate resampling algorithm, and you end up with very grainy map at those lower resolutions. Gdaladdo has several more algorithms, with gauss in many cases the most efficient. However, even with the gauss filter, maps resampled at very low resolutions turn out too hazy. For an example see the 1930 map of the Netherlands: http://mapserver.sara.nl/topo/triang/ If you zoom in to more detailed levels, you'll see the way the image sharpens. For the effect on black-white image choose the TMK-map (1850) with the top-center button. Filters like "unsharp mask" would perform much better in these cases. So I would like to propose two enhancements to gdal: - add additional filters to gdalwarp, gauss and the filters mentioned below. - implement more filters for gdalwarp and gdaladdo, e.g. "unsharp mask", or the "mode" filter asked by Jack below. Perhaps even add the possibility to specify parameters, like in ImageMagick I don't know how difficult is, and whether the gdal devs would find this really an improvement for gdal. I could do this with some ImageMagick or Gimp scripts, but it would be a kludge. As we are going to georeference the complete cadastral and topographical map base of the Netherlands from 1832 to 1994 the next few years (millions of map-scans), this exentsion of gdal would come in very handy. And funding it will really be no problem at all. I would appreciate any kind of comment on this, positive or negative. Regards, Jan Dr. J. Hartmann Department of Geography University of Amsterdam On 03/27/2013 10:13 PM, John Twilley wrote: I'm interested in this feature request as well. Adding the mode resampling algorithm to gdalwarp would be very beneficial to my projects, right up there with being able to access the warp API from Python. Is this at all possible? Should I submit a feature request on Trac, or what? Just let me know! Jack. -- mathuin at gmail dot com On Mon, Feb 18, 2013 at 8:17 AM, Jan Hartmann <mailto:j.l.h.hartm...@uva.nl>> wrote: Hi devs, Would it be possible to add gauss and other interpolations to gdalwarp? At the moment I georeference large scans to 2000*2000 tiles at the most detailed scale, and then create 2000*2000 tiles at resolutions of 2, 4 6 etc times the original scale, using gdaladdo and gauss or other interpolations. It would help immensely if I could do that directly with gdalwarp. Funding would probably no problem. The question is: can and should it be done? Jan ___ gdal-dev mailing list gdal-dev@lists.osgeo.org <mailto:gdal-dev@lists.osgeo.org> http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Ignore geometry_columns table in VRT
Tyler, Hmm. I did the same thing with spatial data in a non-spatial version of SQL Server (via ODBC), using a VRT within a .map file. Not sure if I had any polygons; maybe just points and lines. Must be a bug Best Regards, Brent Fraser On 12/31/2012 3:59 PM, Tyler Mitchell wrote: On 2012-12-31, at 2:03 PM, Brent Fraser wrote: Tyler, So using in the VRT doesn't work when used with Mapserver? Hi Brent, Sort of yes, more or less, but upstreams errors from my OGR driver are somehow stopping it. The main issue is with using a spatial db driver, on a non-spatial db - OGR wants to read the geometry_columns table but it doesn't exist. That's the error that's stopping me. So I'll chalk it up to a bad driver and work on fixing it that way. Thanks! Tyler ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Ignore geometry_columns table in VRT
Tyler, So using in the VRT doesn't work when used with Mapserver? Best Regards, Brent Fraser On 12/31/2012 2:04 PM, Tyler Mitchell wrote: On 2012-12-31, at 12:17 PM, Stephen Woodbridge wrote: Just from code looking, I see that the geometry_columns table is queried from OGRIngresLayer::FetchSRSId(OGRFeatureDefn *poDefn), which is called by OGRIngresTableLayer::ReadTableDefinition(). There's no way (apart changing the code of course!) to disable that query. But I'm not clear why a failure of it would cause MapServer to misbehave if ogrinfo works. I think Mapserver checks this to see if it needs to reproject the data. You might be able to get around it by adding or removing "using srid " clause from the sql. Basically you want to trick mapserver into thinking that it does not need to do a reprojection then it might not make the check that is failing. Mind you this is only my intuition, I have not looked at the code. Thanks guys, I'll give it a whirl. Obviously I need to tweak our Ingres driver so it can work better too :) The weirdest part was having different behaviour between points/lines and the polygons.. but I've found a workaround for now - just used a spatially enabled db instead heh. Tyler ___ 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] OGR XYZ -> GDAL VRT?
Hmmm. You need the GDAL-VRT format to have the same data-source flexibility currently in the OGR-VRT. Could be a useful enhancement with the growing interest in storing rasters in a RDBMs... Best Regards, Brent Fraser On 12/20/2012 10:38 AM, Tyler Mitchell wrote: On 2012-12-20, at 7:34 AM, Brent Fraser wrote: I've also constructed views in the database, basically reformatting the lat/lon columns into WKT format to do: Hi Brent, you gave some good examples, thanks. This has been the approach I took so far, but is not precisely what I'm after. Hopefully I've understood at least some of what you are trying to accomplish. At first I thought you wanted a GDAL VRT that defines a raster, referencing an OGR VRT as the source data and defining a method of producing a grid. Yikes! Yeah, that's the ticket actually :) Think of it as poor-man's raster-in-db solution - a table with X, Y, Z columns - I'm looking for the best way to let GDAL pull that data out and re-assemble it as a raster. I can do it programatically but want to do it on the fly, so MapServer, QGIS, etc. can simply use it. Thanks again, Tyler ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR XYZ -> GDAL VRT?
Tyler, A while back I had a need to access (by Mapserver) points stored in a non-spatial SQL Server database. I was able to use OGR's VRT to access the data, and while the example below shows a datasource connection of ODBC, you may be able to use an Oracle OCI connection instead: ODBC:auser/apass@dsname,Wells_view SELECT * FROM Wells_view WHERE SurfaceLongitude > -117.684 AND SurfaceLongitude < -117.682 AND SurfaceLatitude > 52.499 AND SurfaceLatitude < 52.5 wkbPoint NAD83 y='SurfaceLatitude' x='SurfaceLongitude'/> Actually, I'm a little confused by the GeometryType and GeometryField tags in the above example; you may need to experiment. I've also constructed views in the database, basically reformatting the lat/lon columns into WKT format to do: ODBC:auser/apass@dsname,vWM_UserGraphics SELECT WKT,PropertyID,GeomType,GeomID,LayerName,ClassName,Label FROM vWM_UserGraphics WHERE PropertyID=%PropertyID% ANDGeomType = 0 geomID wkbLineString NAD83 Note the "%PropertyID%" string is a Mapserver runtime substitution done before OGR gets the VRT definition so not very useful in your case (but it might give you a few ideas). Hopefully I've understood at least some of what you are trying to accomplish. At first I thought you wanted a GDAL VRT that defines a raster, referencing an OGR VRT as the source data and defining a method of producing a grid. Yikes! Best Regards, Brent Fraser On 12/20/2012 12:50 AM, Tyler Mitchell wrote: On 2012-12-19, at 10:54 PM, Even Rouault wrote: Le jeudi 20 décembre 2012 02:08:01, vous avez écrit : Hi all, does anyone know a way to put up a GDAL VRT fed by an OGR datasource that has data in XYZ format (regular gridded) . It doesn't have to ultimately be a VRT, just some similar on-the-fly raster format method, so I can share raster data directly from the OGR datasource. I've been digging around and saw wkt raster/postgis raster but it's really a little different than this approach as it doesn't use arrays. I can share more but thought I'd keep it simple first ;-) Tyler, There's a GDAL XYZ driver, but it has very strict requirements : http://gdal.org/frmt_xyz.html Otherwise you can use the OGR CSV driver (make sure the file extension is .csv), wrap it in a OGR VRT, and then ingest that in gdal_grid to produce a raster. But that's not really on-the-fly. Hi Even, Thanks for the tips I'm so close .. thanks to the XYZ driver which I've been using a lot lately, since I can input/output CSV and both GDAL and OGR can read it, but yeah, GDAL can't read it directly from the db source. My challenge is now to make it pull the CSV from the db in real-time via a VRT. (Actually I'm sure I can do it with a SQL statement in a VRT, my db is just not cooperating at the moment, since it's a non-spatial version of the Ingres db and the OGR driver is looking for geometry_columns table and bails out.) But, yeah it's still not on-the-fly from a GDAL raster perspective. I assume some of the other drivers (DODS?) do have a similar mechanism to deal with grids in databases? Isn't there some fancy bash scripting that can be done to masquerade an OGR VRT/script as a virtual file - e.g. it looks like a data file but it's really a realtime script inside it. Which leads me to wonder if there are ways to pipe CSV data into gdalinfo/translate? Hm... thanks again, I'll keep pick away on it when I have time. A generic gdal-ogr CSV bridge driver would be cool! Some for rambling, it's getting late ;) Tyler ___ 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 a shape from another with attributes and values
You need to do a CreateField() on your destination, something like: poFDefn = poLayer->GetLayerDefn(); for( int iField = 0; iField < poFDefn->GetFieldCount(); iField++ ){ OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField ); if( poDstLayer->CreateField( poFieldDefn ) != OGRERR_NONE ){ printf( "Creation of field failed.\n" ); } } Best Regards, Brent Fraser On 12/4/2012 4:53 AM, SIVA RAMA KRISHNA wrote: To All, I am using the following code for creating a shape file from another with the attributes value; i am unable to fetch attribute values from source and change the features from newly created shape files Any Sort of help is greatly appriciated #include "main.h" #include #include #include //OGRFeature* TranslateFeature(OGRFeature* poSrcFeature,OGRFeatureDefn *poFeatureDefn); OGRLayer*poSrcLayer; OGRLayer*poDstLayer; OGRDataSource *poDS; OGRDataSource *poODS=NULL; OGRFeature *poDstFeature; OGRFeature *poSrcFeature; OGRFeatureDefn *poFeatureDefn; int main() { int iField; const char *pszDriverName = "ESRI Shapefile"; OGRSFDriver *poDriver; OGRRegisterAll(); QListfieldid; QListfieldname; const char* str="indiaaa" ; poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName( pszDriverName ); if( poDriver == NULL ) { printf( "%s driver not available.\n", pszDriverName ); exit( 1 ); } poDS = poDriver->Open( "/home/support/Images/IND_adm (3)/IND_adm1.shp" ); poSrcLayer=poDS->GetLayerByName("IND_adm1"); poDriver->DeleteDataSource("firsts1.shp"); poODS = poDriver->CreateDataSource( "firsts1.shp", NULL ); poDstLayer = poODS->CreateLayer( "firsts1",NULL,wkbPolygon); poSrcLayer->ResetReading(); int i; OGRGeometry *poGeometry; // OGRGeometry *bufferGeometry; poDstFeature = OGRFeature::CreateFeature(poDstLayer->GetLayerDefn()); poDstLayer->CreateFeature(poDstFeature); poSrcFeature=poSrcLayer->GetNextFeature(); do{ i=1; OGRFeatureDefn *poFDefn = poSrcLayer->GetLayerDefn(); int iField; for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ ) { qDebug()<<"fieldcount"<GetFieldCount(); printf("\n"); OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField ); if( poFieldDefn->GetType() == OFTInteger ) { printf( "%d,", poSrcFeature->GetFieldAsInteger( iField ) ); } else if( poFieldDefn->GetType() == OFTReal ) { printf( "%.3f,", poSrcFeature->GetFieldAsDouble(iField) ); fieldid.append(poSrcFeature->GetFieldAsDouble( iField )); } else if( poFieldDefn->GetType() == OFTString ) { printf( "%s, Name", poSrcFeature->GetFieldAsString(iField) ); fieldname.append(poSrcFeature->GetFieldAsString(iField)); } else { printf( "%s,name1", poSrcFeature->GetFieldAsString(iField) ); fieldname.append(poSrcFeature->GetFieldAsString(iField)); } } poGeometry=poSrcFeature->GetGeometryRef(); poDstFeature->SetGeometryDirectly(poGeometry); poDstLayer->CreateFeature(poDstFeature); poDstFeature->SetFID(100); // poDstFeature->SetField(i,str); // this is not modifying // i++; // //poDstFeature->SetFrom(poSrcFeature,true); }while( (poSrcFeature = poSrcLayer->GetNextFeature()) != NULL ); poDstLayer->ResetReading(); qDebug()<<"size"<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] Orthorectification iN GDAL
As far as I know, GDAL will orthorectify images supplied with RPC data (e.g. Digital Globe's Worldview/Quickbird, GeoEye's Iknonoe, GeoEye1, etc) but not Spot imagery. You could use OSSIM to orthorectify Spot. Best Regards, Brent Fraser On 10/17/2012 5:22 AM, Luis Lisboa wrote: Greetings I want to orthorectify a few SPOT images and at QGIS mailing list someone indicated me that this feature is available in GDAL. I have checked but I didn't find anything. Can anyone confirm me this? Regards Luis ___ 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 GIS Raster Viewer
Hey Tyler! OSSIM can be a heavy (but useful) raster toolkit, as can OTB. I think there is room in the FOSS4G world for a small GUI app to do raster operations, basically supplying a GUI to gdalwarp, gdal_translate, etc. Kinda like ogr2gui (http://www.ogr2gui.ca) but for raster... Best Regards, Brent Fraser On 10/12/2012 10:37 AM, Tyler Mitchell wrote: I'd recommend the OSSIM platform which as various additional image processing utilities that you will find useful. If you just need desktop GIS viewing functionality, I suggest QGIS which will read all GDAL image formats fine and provide all of what you talk about. Tyler http://trac.osgeo.org/ossim/ http://qgis.org On 2012-10-12, at 5:52 AM, kishore Borra wrote: Hi All, we are implementing support of elevation processing of various raster images in a standalone desktop Viewer using opensource. we are using GDAL library to support performing various operations like hillshade, slope, aspect map & color coded relief. I would like to know your suggestions on what opensource control can be used as a map control so that I can load various type of processed raster images into it. I need to provide the user with the coordinates readout information and standard zoom In, zoom out & zoom By area operations. I should also be able to read the elevation info from the displayed raster and display the same in the status bar. Please suggest. Regards, Kishore. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org <mailto:gdal-dev@lists.osgeo.org> http://lists.osgeo.org/mailman/listinfo/gdal-dev *Tyler Mitchell*** Engineering Director Actian Corporation _tyler.mitch...@actian.com <mailto:tyler.mitch...@actian.com> _ MOBILE 250-303-1831 SKYPE spatialguru www.actian.com <http://www.actian.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] gdaltransform with gcps from file
Jan, Dunno about using the parameter in gdaltransform, but while using gdal_translate, I have used the "--optfile" option and put multiple "-gcp" commands in it. Perhaps it would work with gdaltransform too... Best Regards, Brent Fraser On 6/23/2012 6:23 AM, Jan Hartmann wrote: Silly question, perhaps. I use gdaltransform with gcps given at the command line. According to the documentation, gcps can also be put into the parameter to the program. What is the format of that file? Most of the time I get an error message like "file fomat not supported". With four numbers per line the message is "Ungridded dataset". What am I doing wrong? Jan ___ 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] image outline polygons with keys
All, I need to create image footprints (outlines excluding the nodata areas) for about 300 geotiff images. While GDAL's gdal_polygonize.py may do the job, I prefer GINA's (http://www.gina.alaska.edu/projects/gina-tools) gdal_trace_outline because of its option,s but I'm open to alternatives. I intend to load the polygons into PostGIS for later processing and I need some kind of key value (filename?) on the polygons. Neither gdal_polygonize nor gdal_trace_outline have an option for setting an attribute value (the produced shapefiles have a single integer attribute) to identify the polygon when they create the shapefile. So my current plan calls for scripting: For 300 shapefiles: Use shp2pgsl to Load shapefile into a PostGIS temp table (creates a table with a single row), Use pgsql to do an ALTER table to: add a IMAGENAME column update IMAGENAME for the loaded row copy the row into the final table delete the temp table I suppose I could enhance gdal_polygonize (or gdal_trace_outline) to allow an attribute value, but I'm hoping someone has a simpler, more elegant way using existing software. GeoKettle may be an answer, but it has a steep learning curve... -- Best Regards, Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] JAXAPalsar 1.0 & 1.1 extended metadata reading support
Has there been any progress in reading AVNIR2 data? GDAL's ceos driver doesn't support georeferencing... Best Regards, Brent Fraser On 6/27/2011 6:56 AM, Rodolfo Bonnin wrote: Hello Antonio, I'll read the specs regarding JaxaPalser type, out software is reading those and they look well, but I didn't analyze the format in detail. I've also added georeferencing (Calculating the affine transform based on the 4 corners and assigning a coordinate system (WGS84 instead of ITRF97 at the moment as we are based on 1.8.0) to the JaxaPalsar driver. For the time being, I'm writing new drivers for the PRISM and AVNIR2 sensors, based on the current JaxaPalsar driver, as we are constrained in time, but I think that once I send the code, the diffing could be useful when writing the new unified driver. I'm ready to collaborate if you think another approach could more useful. The idea is to have PRISM and AVNIR2 support ASAP. Regards, Rodolfo Bonnin. 2011/6/26 Antonio Valentino <mailto:antonio.valent...@tiscali.it>>: > Hi Frank, > > Il 26/06/2011 13:42, Frank Warmerdam ha scritto: >> On 11-06-26 06:34 AM, Antonio Valentino wrote: >>> Hi Frank, > > [...] > >>> Looking deeper into the ceos2 code it seems to me that the "recipe" for >>> decoding ALOS-PALSAR products do not match at all specifications I can >>> find at >>> >>> http://www.eorc.jaxa.jp/ALOS/en/doc/fdata/PALSAR_L10_J_ENa.zip >>> >>> Does that recipe refers to another format specification document? >> >> Antonio, >> >> I don't know. > > well, in this case I guess it is not a good idea to modify that recipe > in any manner > >>> Also, regarding the ceos2 driver in general, I didn't find an obvious >>> way to support CEOS products with multiple image files (like multi >>> polarization PALSAR or ALOS PRISM). >> >> Ah, that may well have been part of Phil's rationale for handling >> multi-polarization data with custom drivers. I don't know how hard >> it would be to do so within the existing SAR_CEOS driver. > > OK, thank you for confirming my analysis > >> Possibly you could set up a multi-polarization driver based on the >> low level ceos2 code but with a somewhat different recipe mechanism? > > Yes I could try. > > A new CEOS multi-band driver could be used for ALOS (all three sensors) > and Landsat. > There are other sensors that could be interesting to support? > Where can I find pointers and format spec? > > ciao > > -- > Antonio Valentino > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org <mailto:gdal-dev@lists.osgeo.org> > http://lists.osgeo.org/mailman/listinfo/gdal-dev > -- Ing. Rodolfo Bonnin SUR Emprendimientos Tecnológicos Perú 345 Piso 5to Oficina "B" (C1067AAG) Ciudad de Buenos Aires, Argentina Tel. +54 (11) 4342-2976/84 rodolfobon...@suremptec.com.ar <mailto:rodolfobon...@suremptec.com.ar> www.suremptec.com <http://www.suremptec.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] RE [Polluriel potentiel] Re: ogr ODBC problem [SOLVED]
Steve, That's good news. I did run into the "spaces" problems on Windows (I had a space between the comma and the table name). It caused the connection to fail. Best Regards, Brent Fraser On 5/31/2012 12:14 PM, steve.tout...@inspq.qc.ca wrote: for ogrinfo we were using, as explained here, http://mapserver.org/input/vector/VirtualSpatialData.html#steps-for-display ogrinfo ODBC:User/Pwd@mssql__archives msStatut_station2002_2006 That was working on windows Though, That doesn't work on our linux distro linux sles 11 sp1. It as to be (comma instead of space) ogrinfo ODBC:User/Pwd@mssql__archives,msStatut_station2002_2006 But, For some reasons, this still doesn't work to get the list of tables on our linux distro ogrinfo ODBC:User/Pwd@mssql__archives But it is useless in my case. Also, our virtual data file looks like this on windows ODBC:User/Pwd@mssql__archives dbo.msStatut_station2002_2006 wkbPoint WGS84 That is working on windows For our distro linux We need to add the table name ODBC:User/Pwd@mssql__archives,msStatut_station2002_2006 Perhaps it is not a windows vs linux issue but an OGR version issue. Thanks for your help Steve *Brent Fraser * 2012-05-31 10:49 A steve.tout...@inspq.qc.ca cc gdal-dev@lists.osgeo.org Objet [Polluriel potentiel] Re: [gdal-dev] ogr ODBC problem Steve, Here's a snippet from a MapServer map file I used to access a non-spatial MS SQL server via OGR+VRT+ODBC: TYPE POINT CONNECTIONTYPE OGR CONNECTION " ODBC:username/password@WebMapDSN,v_MyPoints SELECT * FROM v_MyPoints WHERE PropertyID=%PropertyID% *PointID* wkbPoint x='SurfaceLongitude'/> NAD83 " PROCESSING "CLOSE_CONNECTION=DEFER" DATA "v_MyPoints" In my case I was able to use the view's primary key as the FID. Hope this helps... Best Regards, Brent Fraser On 5/31/2012 7:48 AM, _steve.tout...@inspq.qc.ca_ <mailto:steve.tout...@inspq.qc.ca>wrote: Thanks Jeff I got now OGR_ODBC: Table ?s???s!.? has no identified FID column. I found that several users had this problem but found no solution. I don't have write access to this MSSQL server. I'm connecting via ODBC to a non spatial table, but it contains latitude and longitude information. I will use it to define a OGRVRTDataSourceand create geometry from point. Any clue on what I can do? thanks steve *Jeff McKenna **__* <mailto:jmcke...@gatewaygeomatics.com>*@lists.osgeo.org* Envoyé par : _gdal-dev-bounces@lists.osgeo.org_ <mailto:gdal-dev-boun...@lists.osgeo.org> 2012-05-30 17:06 A _gdal-dev@lists.osgeo.org_ <mailto:gdal-dev@lists.osgeo.org> cc Objet Re: [gdal-dev] ogr ODBC problem On 12-05-30 5:09 PM, _steve.tout...@inspq.qc.ca_ <mailto:steve.tout...@inspq.qc.ca>wrote: > > Hi! > I use this command to get the tables from an ODBC connection > ogrinfo ODBC:User/Pwd@DNS > > The connection is succesful but I get this error several times > ERROR 1: No column definitions found for table '?s???s!.???', > layer not usable. > > I used OGR ODBC for several months from a Windows server to a MSSQL server > Now I'm migrating to linux and accessing the same MSSQL SERVER and I get > this error. > > Note that if I connect with isql, I can connect and query the database > without problem. > So the problem really seems to be with OGR > > What can cause this No column definitions found for table error > And why the table name looks like this '?s???s!.???', Hi Steve, I was just debugging an ogrinfo command for MS4W/Oracle, and was reminded of the trick to show more debug info at the commandline: set CPL_DEBUG=on I'm going to keep that one in my back pocket from now on. -jeff -- Jeff McKenna MapServer Consulting and Training Services_ __http://www.gatewaygeomatics.com/_ ___ gdal-dev mailing list_ __gdal-dev@lists.osgeo.org_ <mailto:gdal-dev@lists.osgeo.org>_ __http://lists.osgeo.org/mailman/listinfo/gdal-dev_ ___ gdal-dev mailing list _gdal-dev@lists.osgeo.org_ <mailto: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 [Polluriel potentiel] Re: ogr ODBC problem
Steve, In my case I was getting OGR to get data from a database view, but it could have been a table. I think OGR needs (would like?) to know which column it can use for unque id. I would expect your view (or table) would have a unique key column otherwise relation operations (like joins) can be a problem (the db admin should be able to tell you). Perhaps the ODBC system on Linux doesn't give OGR enough info about the columns to automatically identify the Primary Key column. If you are testing with ogrinfo you could user the "-fid" option after you've identified the Primary Key column? Best Regards, Brent Fraser On 5/31/2012 10:11 AM, steve.tout...@inspq.qc.ca wrote: Thanks Brent for your answer. But view's primary key? I'm missing something here...What do you mean by that? Do you mean a column with a unique id in the view? There is no such column in my view. Is ogr now needs a unique id in a view? Because it was working before I migrate to linux with a recent vertsion of gdal. Also, My problem is before using CONNECTIONTYPE OGR CONNECTION " A simple ogrinfo like this ODBC:username/password@WebMapDSN,v_MyPoints Gives me the error OGR_ODBC: Table ?s???s!.? has no identified FID column. Thanks in advance Steve *Brent Fraser * 2012-05-31 10:49 A steve.tout...@inspq.qc.ca cc gdal-dev@lists.osgeo.org Objet [Polluriel potentiel] Re: [gdal-dev] ogr ODBC problem Steve, Here's a snippet from a MapServer map file I used to access a non-spatial MS SQL server via OGR+VRT+ODBC: TYPE POINT CONNECTIONTYPE OGR CONNECTION " ODBC:username/password@WebMapDSN,v_MyPoints SELECT * FROM v_MyPoints WHERE PropertyID=%PropertyID% *PointID* wkbPoint x='SurfaceLongitude'/> NAD83 " PROCESSING "CLOSE_CONNECTION=DEFER" DATA "v_MyPoints" In my case I was able to use the view's primary key as the FID. Hope this helps... Best Regards, Brent Fraser On 5/31/2012 7:48 AM, _steve.tout...@inspq.qc.ca_ <mailto:steve.tout...@inspq.qc.ca>wrote: Thanks Jeff I got now OGR_ODBC: Table ?s???s!.? has no identified FID column. I found that several users had this problem but found no solution. I don't have write access to this MSSQL server. I'm connecting via ODBC to a non spatial table, but it contains latitude and longitude information. I will use it to define a OGRVRTDataSourceand create geometry from point. Any clue on what I can do? thanks steve *Jeff McKenna **__* <mailto:jmcke...@gatewaygeomatics.com>*@lists.osgeo.org* Envoyé par : _gdal-dev-bounces@lists.osgeo.org_ <mailto:gdal-dev-boun...@lists.osgeo.org> 2012-05-30 17:06 A _gdal-dev@lists.osgeo.org_ <mailto:gdal-dev@lists.osgeo.org> cc Objet Re: [gdal-dev] ogr ODBC problem On 12-05-30 5:09 PM, _steve.tout...@inspq.qc.ca_ <mailto:steve.tout...@inspq.qc.ca>wrote: > > Hi! > I use this command to get the tables from an ODBC connection > ogrinfo ODBC:User/Pwd@DNS > > The connection is succesful but I get this error several times > ERROR 1: No column definitions found for table '?s???s!.???', > layer not usable. > > I used OGR ODBC for several months from a Windows server to a MSSQL server > Now I'm migrating to linux and accessing the same MSSQL SERVER and I get > this error. > > Note that if I connect with isql, I can connect and query the database > without problem. > So the problem really seems to be with OGR > > What can cause this No column definitions found for table error > And why the table name looks like this '?s???s!.???', Hi Steve, I was just debugging an ogrinfo command for MS4W/Oracle, and was reminded of the trick to show more debug info at the commandline: set CPL_DEBUG=on I'm going to keep that one in my back pocket from now on. -jeff -- Jeff McKenna MapServer Consulting and Training Services_ __http://www.gatewaygeomatics.com/_ ___ gdal-dev mailing list_ __gdal-dev@lists.osgeo.org_ <mailto:gdal-dev@lists.osgeo.org>_ __http://lists.osgeo.org/mailman/listinfo/gdal-dev_ ___ gdal-dev mailing list _gdal-dev@lists.osgeo.org_ <mailto: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] ogr ODBC problem
Steve, Here's a snippet from a MapServer map file I used to access a non-spatial MS SQL server via OGR+VRT+ODBC: TYPE POINT CONNECTIONTYPE OGR CONNECTION " ODBC:username/password@WebMapDSN,v_MyPoints SELECT * FROM v_MyPoints WHERE PropertyID=%PropertyID% *PointID* wkbPoint x='SurfaceLongitude'/> NAD83 " PROCESSING "CLOSE_CONNECTION=DEFER" DATA "v_MyPoints" In my case I was able to use the view's primary key as the FID. Hope this helps... Best Regards, Brent Fraser On 5/31/2012 7:48 AM, steve.tout...@inspq.qc.ca wrote: Thanks Jeff I got now OGR_ODBC: Table ?s???s!.? has no identified FID column. I found that several users had this problem but found no solution. I don't have write access to this MSSQL server. I'm connecting via ODBC to a non spatial table, but it contains latitude and longitude information. I will use it to define a OGRVRTDataSourceand create geometry from point. Any clue on what I can do? thanks steve *Jeff McKenna @lists.osgeo.org* Envoyé par : gdal-dev-boun...@lists.osgeo.org 2012-05-30 17:06 A gdal-dev@lists.osgeo.org cc Objet Re: [gdal-dev] ogr ODBC problem On 12-05-30 5:09 PM, steve.tout...@inspq.qc.ca wrote: > > Hi! > I use this command to get the tables from an ODBC connection > ogrinfo ODBC:User/Pwd@DNS > > The connection is succesful but I get this error several times > ERROR 1: No column definitions found for table '?s???s!.???', > layer not usable. > > I used OGR ODBC for several months from a Windows server to a MSSQL server > Now I'm migrating to linux and accessing the same MSSQL SERVER and I get > this error. > > Note that if I connect with isql, I can connect and query the database > without problem. > So the problem really seems to be with OGR > > What can cause this No column definitions found for table error > And why the table name looks like this '?s???s!.???', Hi Steve, I was just debugging an ogrinfo command for MS4W/Oracle, and was reminded of the trick to show more debug info at the commandline: set CPL_DEBUG=on I'm going to keep that one in my back pocket from now on. -jeff -- Jeff McKenna MapServer 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 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] how to convert UInt16 Tiff
Another alternative is GINA's (http://www.gina.alaska.edu/projects/gina-tools/) gdal_contrast_stretch: Usage: gdal_contrast_stretch.exe { { -linear-stretch } | { -percentile-range } | { -histeq } } [ -ndv ] Input must be either 8-bit or 16-bit. Output is 8-bit. I think the input can be any GDAL supported raster type. Best Regards, Brent Fraser On 4/24/2012 2:53 AM, jr.morre...@enoreth.net wrote: On Tue, 24 Apr 2012 09:38:44 +0200, Even Rouault wrote: thanks for the reply, I did try with NBITS=8 and 12 with osgeo4w's gdal and tamas' nightly build but the output results are incorrect Here is a sample file (93Mo) : http://dl.free.fr/iFT4hH3Yj I'll try to rebuild gdal on my linux box with the internal libs I haven't looked at your result image, but you cannot just specify NBITS=8 or 12, if the range of the values in the original file is full 16 bits. You need to rescale the values, otherwise they will get clamped to [0, 255] or [0, 4096]. So you can look for the min/max reported by gdalinfo -mm uncompressed_original.tif and then try : gdal_translate -co "COMPRESS=JPEG" uncompressed_original.tif compressed.tif -scale min max 0 255 -ot Byte or gdal_translate -co "COMPRESS=JPEG" uncompressed_original.tif compressed.tif -scale min max 0 4095 -co NBITS=12 Thanks, I didn't knew that ! The man says that you can omit the input min max scale so gdal computes it for each source, however I can't find how to write that correctly ___ 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] Contrast, Brightness and Gamma
I like to use gdal_contrast_stretch from http://www.gina.alaska.edu/projects/gina-tools/ Usage: gdal_contrast_stretch { { -linear-stretch } | { -percentile-range } | { -histeq } } [ -ndv ] Input must be either 8-bit or 16-bit. Output is 8-bit. For example: gdal_contrast_stretch -percentile-range 0.02 0.98 input16bit.tif output8bit.tif Best Regards, Brent Fraser On 3/20/2012 9:55 AM, Saâd HESSANE wrote: Hi list, I use gdal_translate to convert 16bits images to 8bits images. The gdal_translate have the -scale argument to specify the convertion range (source->destination). If I change the destination range (dst_min and dst_max), I can do manually a brightness correction (also a contrast correction). But is there any way to do that with a Photoshop like method? For example set a brightness to -50, a contrast to +8 and let gdal_translate do the job ? Also can I apply a gamma correction? I saw the VRT format and there's a way to apply a LUT. But I don't think a LUT of 65536 values is the best way to apply the correction. Also the VRT have the scale ratio and scale offset elements. I think it can be useful but I don't understand how use it... So simply put, my question : can I do a contrast, brightness or gamma correction to a raster with gdal? Thank you ! PS : I apologize for my broken English! ___ 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: PDF selective layer sample wanted
I need more coffee. Thanks! Best Regards, Brent Fraser On 3/5/2012 8:22 AM, Rahkonen Jukka wrote: Hi, I did mean that with "contains vector graphics" and "zoom in and enjoy". -Jukka- Brent Fraser wrote: Cool. Now if we could get it writing vectors... Best Regards, Brent Fraser On 3/5/2012 8:10 AM, Jukka Rahkonen wrote: Brent Fraser geoanalytic.com> writes: I thought the PDF driver was read-only. Does the changeset include the creation of a geospatial pdf? Read thread "Trouble with georeferenced PDF" from couple of days ago. Here is also a new sample which is created with Mapserver from trunk and it contains vector graphics. Metadata is in OGC_BP format and you'll need TerraGo toolbar for seeing the coordinates. Zoom in and enjoy. http://latuviitta.org/documents/Lahti_geospat.pdf -Jukka- ___ 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: PDF selective layer sample wanted
Cool. Now if we could get it writing vectors... Best Regards, Brent Fraser On 3/5/2012 8:10 AM, Jukka Rahkonen wrote: Brent Fraser geoanalytic.com> writes: I thought the PDF driver was read-only. Does the changeset include the creation of a geospatial pdf? Read thread "Trouble with georeferenced PDF" from couple of days ago. Here is also a new sample which is created with Mapserver from trunk and it contains vector graphics. Metadata is in OGC_BP format and you'll need TerraGo toolbar for seeing the coordinates. Zoom in and enjoy. http://latuviitta.org/documents/Lahti_geospat.pdf -Jukka- ___ 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: PDF selective layer sample wanted
I thought the PDF driver was read-only. Does the changeset include the creation of a geospatial pdf? Best Regards, Brent Fraser On 3/5/2012 7:59 AM, Jukka Rahkonen wrote: Brent Fraser geoanalytic.com> writes: If you're looking for a geospatial pdf with layers, the Canadian government publishes some of their topographic maps as pdf. Here's a link to one: ftp://ftp2.cits.rncan.gc.ca/pub/cantopo/50k_geopdf/085/ b/cantopo_085b14_geopdf.zip Thanks for the answer. I do already have geospatial pdf files with layers. What I am after is geospatial pdf files with layers and created with GDAL. I can now make pdf files with layers with SkyJUMP GIS and geospatial pdf files with GDAL but I got hungry and want them both together. -Jukka- ___ 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] PDF selective layer sample wanted
If you're looking for a geospatial pdf with layers, the Canadian government publishes some of their topographic maps as pdf. Here's a link to one: ftp://ftp2.cits.rncan.gc.ca/pub/cantopo/50k_geopdf/085/b/cantopo_085b14_geopdf.zip Best Regards, Brent Fraser On 3/5/2012 1:29 AM, Jukka Rahkonen wrote: Hi, This changeset looks interesting http://trac.osgeo.org/gdal/changeset/24069 Is there any such PDF sample file available? -Jukka Rahkonen- ___ 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] org2ogr, KML, and srs definition
I have some shapefiles I thought I would look at in Google Earth, so I used ogr2ogr (v1.8.1) to convert them to KML: ogr2ogr -f "KML" Lines_OGR.kml Lines.shp but they appeared about 100 meters Northeast of my expectation, leading me to believe there was a Datum problem. ogrinfo reported (since there is a .prj file created by ArcMap): Layer SRS WKT: PROJCS["UTM_Zone_35_Northern_Hemisphere", GEOGCS["GCS_Geographic Coordinate System", DATUM["EUROPEAN_DATUM_1950", SPHEROID["International_1909_Hayford_Intl_1924",6378388,297.000284015]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",27], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",50], PARAMETER["false_northing",0], UNIT["Meter",1]] note the ED50, which is correct for my shapefiles. After a few attempts experimenting, it was only when I used -s_srs AND -t_srs on the command line the features were shown where I expected them: ogr2ogr -s_srs "EPSG:23035" -t_srs "EPSG:4326"... The part that bothered me was there was no message from ogr2ogr warning me about source/target SRS issues. It seems counter-intuitive since ogrinfo was able to intreprete the .prj file of the source data, and the only (?) allowable SRS of KML is EPSG:4326. It this a known issue (other than a related http://trac.osgeo.org/gdal/ticket/2271)? -- Best Regards, Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GeoPDF translation
Michael, GDAL's PDF format page says you can use the GDAL_PDF_DPI config option... Best Regards, Brent Fraser On 12/13/2011 11:58 AM, Smith, Michael wrote: I am trying to convert the historic USGS topo maps (GeoPDFs) into GeoTIFFs with GDAL 1.8 . A very straightforward thing to do is gdal_translate in.pdf out.tiff This works fine but the TIFF looks as if the PDF was exported at 150dpi, which is a pretty crappy resolution. Using something else (like PDFCreator and printing to a TIFF), the USGS GeoPDFs can be exported at higher resolutions, but then of course you lose the georeferencing. So I have GDAL which keeps the georeferencing but only seems to like 150dpi, and then a zillion other solutions which allow 300 dpi or 600 dpi, but have no georeferencing. I did try gdalwarp -tr 5 5 for example (which would roughly be equivalent to one of these maps at 300dpi), but it is clear that it is still rendering the PDF at 150dpi and then generating a 5-meter TIFF - looks no different really than the default output with gdal_translate. Is there some argument that would specify the dpi at which the GeoPDF were rendered before conversion to GeoTIFF? I don't see anything like that in the docs or the list archives. BTW if you haven't checked out the historic topos they are great, going back to the 1880s for Maine. What fun! My ultimate goal is to have a seamless mosaic WMS of the old topos for Maine. * Michael Smith State GIS Manager, Maine Office of GIS, Maine OIT Board Member, Maine GeoLibrary Board Member, Maine GIS Users Group State Rep, National States Geographic Information Council State House Station 174 264 Civic Center Drive Augusta, ME 04333-0174 (207) 215-5530 69 o 47' 49.5"W 44 o 20' 54.5"N ___ 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] how to get gis layer information from GeoPDF (e.g. retrieve building polygon coordinates)
I wonder if it would be possible to leverage pstoedit (http://www.pstoedit.net/)? I've used it with Ghostview to extract PDF vectors to DXF (with no georeferencing of course) with varying success. Best Regards, Brent Fraser On 12/2/2011 3:37 PM, Even Rouault wrote: Le vendredi 02 décembre 2011 23:19:10, Jared Rubin a écrit : Even, Would it be possible to do an ogr driver to get the vector features (I guess similar to shape files? Or maybe update the pdf driver to dump each indivual vector feature into gdal metadata? A OGR driver would perhaps be possible. I haven't studied the question however. And I presume it would represent a significant effort as it would involve parsing and interpretating the PDF vector stream. Currently, all the hard work of rasterizing the PDF content is delegated to the poppler library. ___ 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: Problem joining a shapefile and MDB
A response from the mapserver list: http://forums.devarticles.com/microsoft-sql-server-5/data-source-name-not-found-and-no-default-driver-specified-8346.html When you use a 64bits Windows, ODBC manager at Control Panel will only create 64bits connection. I had problems with that in the past and had to force 32bits connections using 32bits ODBC manager at C:\Windows\SysWOW64\odbcad32.exe Best Regards, Brent Fraser On 11/18/2011 2:00 PM, boesiii wrote: Even, What is the correct procedure to upgrade a MS4W installation? I downloaded release-1310-gdal-1-8-1-mapserver-6-0-1.zip and I copied OGRINFO.EXE to the current location. After I did this OGRINFO still crashed. Do I need to copy other files? Thanks, Ed -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/Problem-joining-a-shapefile-and-MDB-tp7009404p7009680.html Sent from the GDAL - Dev mailing list archive at Nabble.com. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] 2D KML to 3D KML
Hi all, I've got some 2D KML generated by gdal_polygonize I'd like to turn into 3D KML. Do you have any ogr2ogr techniques to do this? A snip of my 2D KML: -- 1 <LineStyle><color>ffff</color></LineStyle><PolyStyle><fill>0</fill></PolyStyle> 1 106.2125,-6.0 106.2125,-6.01 106.225,-6.01 106.225,-6.02 106.2375,-6.02 106.2375,-6.03 106.25,-6.03 106.25,-6.0 106.2125,-6.0 1 <LineStyle><color>ffff</color></LineStyle><PolyStyle><fill>1</fill><color>ffff</color></PolyStyle> 1 106.2125,-6.0,101 106.2125,-6.01,101 106.225,-6.01,101 106.225,-6.02,101 106.2375,-6.02,101 106.2375,-6.03,101 106.25,-6.03,101 106.25,-6.0,101 106.2125,-6.0,101 -- I've added a "z" to the coordinates (100 plus the "name" value, to boost it above the terrain), plus a little bit of PolyStyle info to set the fill color. The long term solution may be to enhance gdal_polygonize to output 3D vectors, but right now I'm just experimenting with visualizing KML in Google Earth -- Best Regards, Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Percent stretch 16 to 8-bit geotiff conversion
Have a look at http://www.gina.alaska.edu/projects/gina-tools/ Best Regards, Brent Fraser On 9/7/2011 2:22 PM, David Shean wrote: Forgive me if this is not the proper forum (no gdal-users list?) or if this question has been answered in the past - I've been using GDAL for a while now, but just joined the mailing list and I'm relatively new to python. I have large (>2GB uncompressed) UInt16 geotiffs that I need to convert to 8-bit geotiffs. I'd like to use a percent range (e.g. 0.5-99.5%) to determine the new min/max for the linear stretch in gdal_translate. I haven't found any simple way of doing this with the existing GDAL utilities, although it looks like gdalenhance was at least initially designed to do something along these lines. If there is no simple way to do this with existing tools, is there any sample code floating around that might help? If not, I know it will be possible with the histogram/stats functions in gdal, numpy, or PIL - I just want to make sure that I'm not reinventing a slower, clunkier wheel. Thanks. -David --- David Shean Polar Science Center Applied Physics Lab University of Washington ___ 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] adding georeferencing information
Jan, I use: listgeo mygeotiff.tif > mygeotiff.gtt geotifcp -g mygeotiff.gtt nogeotags.tif mynewgeotiff.tif FYI, I've been experimenting with GDAL's vrt format to sharpen my tiffs: : : Red NRT_123.tif 1 BlockXSize="2890" BlockYSize="1" /> 3 -0.111 -0.111 -0.111 -0.111 2 -0.111 -0.111 -0.111 -0.111 : : My kernel tends to brighten the image a little (note the "2" in the kernel) as well as sharpen. To try it on your file, create a VRT file with gdal_translate -of VRT, then use a text editor to change the SimpleSource tags to KernelFilteredSource and add the Kernel. You can then create a sharpened tif by using the VRT as input to gdal_translate, os simply open the VRT in Quantum GIS to view the results. Best Regards, Brent Fraser On 9/5/2011 8:46 AM, Jan Hartmann wrote: Hi all, Is there a quick way to add back georeferencing information to a geotif file that has been processed by Imagemagick? I have lots of georeferenced images that I want to sharpen up a bit, but ImageMagick throws the georeference headers away. Is there a way to transport them from the original file to the derived one? Jan ___ 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: Wildfire mapping
Mark, There are a few Open Source projects you could look into 3D visualization environment: VTP WorldWind OssimPlanet Orthorectification of Video: OSSIM OMAR? Image Processing / Feature Extraction OTB GRASS And I'm sure there are others. The system you describe would be a lot of work, but it does sound interesting! Best Regards, Brent Fraser On 8/15/2011 10:35 PM, Mark Zaller (AerialFireTech) wrote: Hello, I'm pursuing a project to automatically map wildfires, and am looking for people who could contribute to this. I volunteer as an Air Attack pilot over wildfires, and have developed an in cockpit Augmented Reality system that allows us to work at night and over smoke. It combines 3D visualization of GIS (Geographic Information System) data along with live Infrared. I'd like to consider porting to Open Source/Linux and include auto-mapping of spotfires. Here is a proposal presentation about what NASA might fund: http://www.slideshare.net/mszaller/wildfire-ir-and-mapping Please write if you have experience with how GRASS, OSGeo, GDAL, or point me to where can find Open Source people who could help. thanks, -Mark Zaller - TriGeo <http://www.trigeocorp.com/> - (408) 623-4303 , www.AerialFireTech.org <http://www.aerialfiretech.org/> Linkedin Profile <http://www.linkedin.com/in/markzaller> - Business Development, Product Management, Hi-tech Hardware & Software Team Leader ___ 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: How to use gdalinfo for HDF- Windows
Luisa, That sounds like a bug. A possible work-around might be to set your current working directory to D and don't put the drive letter in the filepath: cd D: gdalinfo \data\GL.hdf But this may not work with your desire to use a GRASS script... Best Regards, Brent Fraser On 8/10/2011 2:55 AM, Luisa Peña wrote: Hey Brent I'm using ESA's GLOBCOVER data.In this case I'm trying to set a GRASS script to run gdal directly. I do: gdalinfo C:\data\GL.hdf and I retrieve all basic information plus HDF's layers Driver: HDF4/Hierarchical Data Format Release 4 Files: c:\data\GL.hdf Size is 512, 512 Coordinate System is `' Metadata: HDFEOSVersion=HDFEOS_V2.12 Subdatasets: SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"c:\data\GL.hdf":POSTEL:MEAN_1 SUBDATASET_1_DESC=[1800x1800] MEAN_1 POSTEL (16-bit integer) SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"c:\data\GL.hdf":POSTEL:MEAN_2 SUBDATASET_2_DESC=[1800x1800] MEAN_2 POSTEL (16-bit integer) SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"c:\data\GL.hdf":POSTEL:MEAN_3 SUBDATASET_3_DESC=[1800x1800] MEAN_3 POSTEL (16-bit integer) SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"c:\data\GL.hdf":POSTEL:MEAN_4 SUBDATASET_4_DESC=[1800x1800] MEAN_4 POSTEL (16-bit integer) SUBDATASET_5_NAME=HDF4_EOS:EOS_GRID:"c:\data\GL.hdf":POSTEL:MEAN_5 (...) then: gdalinfo HDF4_EOS:EOS_GRID:"c:\data\GL.hdf":POSTEL:MEAN_1 I get: ERROR 4: 'HDF4_EOS:EOS_GRID:c:\data\GL.hdf:POSTEL:MEAN_1' does not exist in fyle system, and is not recognised as a supported dataset name. gdalinfo failed- unable to open 'HDF4_EOS:EOS_GRID:c:\data\GL.hdf:POSTEL:MEAN_1'' open GDAL Datasets: 1 N DriverIsNULL 512x512x0 But if I do: gdalinfo HDF4_EOS:EOS_GRID:":\data\GL.hdf":POSTEL:MEAN_1 it works (without the C from the driver) The problem is that if I'm trying to to use data in D it does not work. Any suggestions? 2011/8/9 Brent Fraser <mailto:bfra...@geoanalytic.com>> Luisa, I use GDAL to access MODIS data in HDF4 on Windows like this: gdalinfo HDF4_EOS:EOS_SWATH:"G:\Projects\temp\image3\MOD02HKM.A2011186.1640.005.NRT.hdf":MODIS_SWATH_Type_L1B:EV_500_RefSB So I'd suggest using double quotes around your path. Best Regards, Brent Fraser On 8/9/2011 8:24 AM, Luisa Peña wrote: Greetings This seems odd but it does not work with HDF4_EOS:EOS_GRID:C:/filepath only HDF4_EOS:EOS_GRID:filepath And if my MERIS data is not placed in my main hard disk (C:) it does not work. Anyone with some experience on this in WINDOWS? Thanks Luisa ___ gdal-dev mailing list gdal-dev@lists.osgeo.org <mailto: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: How to use gdalinfo for HDF- Windows
Luisa, I use GDAL to access MODIS data in HDF4 on Windows like this: gdalinfo HDF4_EOS:EOS_SWATH:"G:\Projects\temp\image3\MOD02HKM.A2011186.1640.005.NRT.hdf":MODIS_SWATH_Type_L1B:EV_500_RefSB So I'd suggest using double quotes around your path. Best Regards, Brent Fraser On 8/9/2011 8:24 AM, Luisa Peña wrote: Greetings This seems odd but it does not work with HDF4_EOS:EOS_GRID:C:/filepath only HDF4_EOS:EOS_GRID:filepath And if my MERIS data is not placed in my main hard disk (C:) it does not work. Anyone with some experience on this in WINDOWS? Thanks Luisa ___ 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] gdalbuildvrt does not support heterogenous band characteristics
Matt, Likely Quickbird_06m10.tif has more or less bands than the previous files, or band bit depth is different. For example, don't try to include a 1 band 16-bit (e.g. gray scale) tiff in a VRT of 3 band, 8 bits each (e.g. RGB color) tiffs. What does gdalinfo say about this file (and other files in the VRT)? Best Regards, Brent Fraser On 6/14/2011 1:13 PM, Matt Wilkie wrote: Hi Folks, what does this error mean? gdalbuildvrt -input_file_list img-list.txt mosaic.vrt 0...10Warning 6: gdalbuildvrt does not support heterogenous band characteristics . Skipping Quickbird_06m10.tif ...etc. thanks! ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Clipping multi-polygons shapefiles with Geometry.Intersection . Any help to fix out my method?
It may help to supply a small test dataset that causes the problem. Best Regards, Brent Fraser On 5/31/2011 1:28 AM, hajer wrote: No idea about this issue?? Thanks for help ! -- View this message in context: http://osgeo-org.1803224.n2.nabble.com/gdal-dev-Clipping-multi-polygons-shapefiles-with-Geometry-Intersection-Any-help-to-fix-out-my-method-tp6410874p6421682.html Sent from the GDAL - Dev mailing list archive at Nabble.com. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: [mapserver-users] OGRinfo works, but a MapServer test mode=map DOES NOT
Eduardo, What happens when you use ogrinfo to access it via ODBC: ogrinfo ODBC:my_user/my_passwd@my_dsn or if you put the connection in a vrt file and use ogrinfo on it? Best Regards, Brent Fraser On 5/26/2011 4:26 AM, Eduardo Kanegae wrote: Hi all, I´m running Windows Server 2008 R2 Enterprise SP1 64bits using MS4W 3.0.1 with MapServer CGI 5.6.6 grabed from MapTools.org My spatial data is on SQL Server 2008 + MsSqlSpatial 0.1.1 Lib and I´ve set an ODBC DSN connection using SQL Server Native Driver 10 ( using Control Panel -> ODBC Admin -> System DSN and also C:\Windows\SysWOW64\odbcad32.exe tool ) When I run ogrinfo using: $ ogrinfo -so "MSSQL:dsn=my_dsn;server=(local);database=my_database;tables=dbo.my_table;uid=my_user;pwd=my_passwd" my_table it works and dumps info like: ... INFO: Open of `MSSQL:dsn=nucoffee;server=(local);database=nucoffee;tables=dbo.fazendas2;uid=db_nucoffee;pwd=@&NuCoFfEe&@' using driver `MSSQLSpatial' successful. ... and so on! BUT, when trying a MapServer CGI mode map request like: http://my_server/cgi-bin/mapserv.exe?map=e:/my_mapfile.map&mode=map&layers=all&map_size=800+600 I only got the error: msDrawMap(): Image handling error. Failed to draw layer named 'pol_my_layer'. msOGRFileOpen(): OGR error. Open failed for OGR connection in layer `pol_my_layer'. Unable to initialize ODBC connection to DSN for my_user/my_passwd@my_dsn, [Microsoft][ODBC Driver Manager] Data source name not found and no default driver specified Can anybody give me a clue on that? My LAYER mapfile part is defined as : ** LAYER NAME"pol_my_layer" TYPEPOLYGON CONNECTIONTYPE OGR CONNECTION " ODBC:my_user/my_passwd@my_dsn dbo.vw_pol_my_spatial_view WGS84 wkbUnknown " DATA "pol_my_layer" FILTER "329" FILTERITEM "farm_id" PROCESSING "CLOSE_CONNECTION=DEFER" END ... END ** best regards, -- Eduardo Patto Kanegae ___ mapserver-users mailing list mapserver-us...@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/mapserver-users ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] validate/clean geometry with OGR?
Elijah, I've used GEOS Buffer function available via the OGR API (http://www.gdal.org/ogr/classOGRGeometry.html) to clean polygons: OGRGeometry *poTempGeometry = poGeometry->Buffer(0.0); // Let Geos fix the problem. Best Regards, Brent Fraser On 4/27/2011 8:36 AM, Elijah Robison wrote: Hey devs, is anyone aware of an OGR approach to validate/clean geometries, for instance, as they are being converted from SHP to PostGRESql? It's not uncommon to have a handful of invalid geometries (usually self-intersecting polygons) in a parcel dataset, and it's impractical to correct them manually. In most situations I use the "Clean Geometries" function in ArcGIS, which does fix the problem. However this one functionality is my main reason to retain Arc in my primary tool chain. I've never found a scripting approach by trolling the search engines so I thought I'd ask here. Thanks, Elijah ___ 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] PDF-->JPEG-->KML. Confused about association between KML bounding box and PDF's neatline and corner coordinates.
Hmmm. My PDF doesn't seem to have a neatline, only a MediaBox and two BBoxes... Brent On 10/20/2010 9:04 AM, Even Rouault wrote: You don't need to convince anyone to add new functionnality. It already exists ;-) See the -cutline option of gdalwarp. Basically you use the NEATLINE reported by gdalinfo as the CUTLINE (that was the sole purpose of reporting the NEATLINE by the way !) The easiest way is to make a simple CSV like this : foo,WKT bla,"POLYGON(())" and you paste the content of NEATLINE as the value of the WKT. and do a gdalwarp without reprojecting to crop (use -crop_to_cutline) to the neatline. Then you can do a 2nd gdalwarp to reproject to the target SRS. (if you want to do just one gdalwarp, you'll have to do ogr2ogr to reproject the cutline.csv into the target SRS) See http://gdal.org/gdalwarp.html and/or http://gdal.org/structGDALWarpOptions.html#0ed77f9917bb96c7a9aabd73d4d06e08 Boris, Use gdalwarp to convert your jpeg to Geographic coordinates (I think it needs to be converted for KML anyways): gdalwarp -t_srs EPSG:4326 in.jpg out.jpg then use gdalinfo to get the corner coordinates. Currently the entire page of the PDF is rendered into the image. Perhaps we can convince Even to add an option to convert only the mapped area (or mask everything else?) Best Regards, Brent Fraser On 10/20/2010 8:07 AM, Boris Dev wrote: My objective is to make a KML ground overlay with a geospatial PDF map. Using GDAL's translate utility I have converted the PDF to a JPEG. The next step and problem for me is to write the text of the KML ground overlay, which requires a text of bounding box information like this: 37.91904192681665 37.46543388598137 15.35832653742206 14.60128369746704 -0.1556640799496235 As GDAL's output, I saw that the entire PDF with legend, title, etc and not just the map image was converted. As GDAL's output, I also saw 2 parts of data that seem they are useful to match a bounding box to the map within the jpeg, but I cant figure out the puzzle of how to make this work? Can I get the bounding box of the entire PDF image with this information below? Corner Coordinates: Upper Left ( 666937.758, 2024687.983) Lower Left ( 666991.811, 2018338.153) Upper Right ( 674511.276, 2024754.795) Lower Right ( 674565.328, 2018404.965) Center ( 670751.543, 2021546.474) POLYGON ((672676.488065659650601 2019417.955749646294862,672636.091432046378031 2024030.267778463661671,667347.405745737953112 2023983.612156898481771,667385.539170471020043 2019370.524140953086317,672676.488065659650601 2019417.955749646294862,672676.488065659650601 2019417.955749646294862,672676.488065659650601 2019417.955749646294862,672676.488065659650601 2019417.955749646294862)) I realize the coordinate above need to be converted to lat/long to be a KML ground overlay onto Google Maps, which I think I can handle. But I am unsure if I need to extract the map image from the jpeg using the NEATLINE? How are the Corner Coordinates translated into a bounding box? Any suggestions? ___ 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] PDF-->JPEG-->KML. Confused about association between KML bounding box and PDF's neatline and corner coordinates.
Boris, Use gdalwarp to convert your jpeg to Geographic coordinates (I think it needs to be converted for KML anyways): gdalwarp -t_srs EPSG:4326 in.jpg out.jpg then use gdalinfo to get the corner coordinates. Currently the entire page of the PDF is rendered into the image. Perhaps we can convince Even to add an option to convert only the mapped area (or mask everything else?) Best Regards, Brent Fraser On 10/20/2010 8:07 AM, Boris Dev wrote: My objective is to make a KML ground overlay with a geospatial PDF map. Using GDAL's translate utility I have converted the PDF to a JPEG. The next step and problem for me is to write the text of the KML ground overlay, which requires a text of bounding box information like this: 37.91904192681665 37.46543388598137 15.35832653742206 14.60128369746704 -0.1556640799496235 As GDAL's output, I saw that the entire PDF with legend, title, etc and not just the map image was converted. As GDAL's output, I also saw 2 parts of data that seem they are useful to match a bounding box to the map within the jpeg, but I cant figure out the puzzle of how to make this work? Can I get the bounding box of the entire PDF image with this information below? Corner Coordinates: Upper Left ( 666937.758, 2024687.983) Lower Left ( 666991.811, 2018338.153) Upper Right ( 674511.276, 2024754.795) Lower Right ( 674565.328, 2018404.965) Center ( 670751.543, 2021546.474) POLYGON ((672676.488065659650601 2019417.955749646294862,672636.091432046378031 2024030.267778463661671,667347.405745737953112 2023983.612156898481771,667385.539170471020043 2019370.524140953086317,672676.488065659650601 2019417.955749646294862,672676.488065659650601 2019417.955749646294862,672676.488065659650601 2019417.955749646294862,672676.488065659650601 2019417.955749646294862)) I realize the coordinate above need to be converted to lat/long to be a KML ground overlay onto Google Maps, which I think I can handle. But I am unsure if I need to extract the map image from the jpeg using the NEATLINE? How are the Corner Coordinates translated into a bounding box? Any suggestions? ___ 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] Re: Testing Geospatial PDF support
lorInterp=Red Band 2 Block=7179x1 Type=Byte, ColorInterp=Green Band 3 Block=7179x1 Type=Byte, ColorInterp=Blue Doing the math on a 1:75k map, at 150 dpi the pixel size should be about 12.7 meters. Brent On 10/19/2010 12:45 PM, Even Rouault wrote: Brent, I'm not sure what you find is wrong. What does gdalinfo returns ? The PDF driver will reproject the coordinates of the GPTS array from geographic into UTM 35N and try to build a geotransform matrix from that. The matrix may have rotating terms, so you likely need to use gdalwarp afterwards to make it look "straight". However I'm not 100% sure this is the correct way to interpret the the Adobe style georeferencing, but I didn't find a better one up to know. So if you or someone else have any clue... Most geospatial PDFs I've found in the wild (such as the one of cantopo) seem to use the OGC Best Practice encoding. Even Le mardi 19 octobre 2010 20:11:13, Brent Fraser a écrit : Even, I've tested a couple of Geospatial PDFs with the code that's in GDAL trunk. The Canadian topos work well (e.g. ftp://ftp2.cits.rncan.gc.ca/pub/cantopo/50k_geopdf/085/b/cantopo_085b14_geo pdf.zip). They were produced using ESRI ArcMap 9.2.6.1500. I get a properly geo-referenced GeoTiff after running gdal_translate. I have another PDF produced by ArcMap10.0.0.2414. It has a main map and a keymap. The strange thing is that even though the coordinate system is UTM Zone 35 , it has a transformation matrix with geographic coordinates given. Here's a portion of the "dump" file: Item[4] : VP Type = array Item[0]: Type = dictionary Item[0] : Name = ÿþI (string) Item[1] : Measure Type = dictionary Item[0] : Subtype = GEO (name) Item[1] : LPTS Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : GPTS Type = array Item[0]: 38.815000 (real) Item[1]: 23.644900 (real) Item[2]: 43.662600 (real) Item[3]: 23.387000 (real) Item[4]: 43.615000 (real) Item[5]: 31.889000 (real) Item[6]: 38.774900 (real) Item[7]: 31.540700 (real) Item[3] : GCS Type = dictionary Item[0] : WKT = PROJCS["WGS_1984_UTM_Zone_35N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPH EROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["De gree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["Fal se_Easting",50.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Me ridian",27.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origi n",0.0],UNIT["Meter",1.0]] (string) Item[1] : Type = PROJCS (name) Item[4] : Type = Measure (name) Item[5] : Bounds Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : Type = Viewport (name) Item[3] : BBox Type = array Item[0]: 3168 (int) Item[1]: 2433 (int) Item[2]: 3557 (int) Item[3]: 2127 (int) Item[1]: Type = dictionary Item[0] : Name = ÿþA (string) Item[1] : Measure Type = dictionary Item[0] : Subtype = GEO (name) Item[1] : LPTS Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : GPTS Type = array Item[0]: 40.190600 (real) Item[1]: 26.279300 (real) Item[2]: 40.767200 (real) Item[3]: 26.273100 (real) Item[4]: 40.769300 (real) Item[5]: 27.212000 (real) Item[6]: 40.192600 (real) Item[7]: 27.210200 (real) Item[3] : GCS Type = dictionary Item[0] : WKT = PROJCS["WGS_1984_UTM_Zone_35N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPH EROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["De gree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["Fal se_Easting",50.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Me ridian",27.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origi n",0.0],UNIT["Meter",1.0]] (string) Item[1] : Type = PROJCS (name) Item[4] : Type = Measure (name) Item[5] : Bounds Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : Type = Viewport (name) Item[3] : BBox Type = array Item[0]: 57 (int) Item[1]: 2465 (int) Item[2]: 3052 (int) Item[3]: 46 (int) The resulting GeoTiff has a scaling problem, I suspect related to the transformation matrix. Any thoughts on that? Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Testing Geospatial PDF support
Even, I've tested a couple of Geospatial PDFs with the code that's in GDAL trunk. The Canadian topos work well (e.g. ftp://ftp2.cits.rncan.gc.ca/pub/cantopo/50k_geopdf/085/b/cantopo_085b14_geopdf.zip). They were produced using ESRI ArcMap 9.2.6.1500. I get a properly geo-referenced GeoTiff after running gdal_translate. I have another PDF produced by ArcMap10.0.0.2414. It has a main map and a keymap. The strange thing is that even though the coordinate system is UTM Zone 35 , it has a transformation matrix with geographic coordinates given. Here's a portion of the "dump" file: Item[4] : VP Type = array Item[0]: Type = dictionary Item[0] : Name = ÿþI (string) Item[1] : Measure Type = dictionary Item[0] : Subtype = GEO (name) Item[1] : LPTS Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : GPTS Type = array Item[0]: 38.815000 (real) Item[1]: 23.644900 (real) Item[2]: 43.662600 (real) Item[3]: 23.387000 (real) Item[4]: 43.615000 (real) Item[5]: 31.889000 (real) Item[6]: 38.774900 (real) Item[7]: 31.540700 (real) Item[3] : GCS Type = dictionary Item[0] : WKT = PROJCS["WGS_1984_UTM_Zone_35N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",50.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",27.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] (string) Item[1] : Type = PROJCS (name) Item[4] : Type = Measure (name) Item[5] : Bounds Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : Type = Viewport (name) Item[3] : BBox Type = array Item[0]: 3168 (int) Item[1]: 2433 (int) Item[2]: 3557 (int) Item[3]: 2127 (int) Item[1]: Type = dictionary Item[0] : Name = ÿþA (string) Item[1] : Measure Type = dictionary Item[0] : Subtype = GEO (name) Item[1] : LPTS Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : GPTS Type = array Item[0]: 40.190600 (real) Item[1]: 26.279300 (real) Item[2]: 40.767200 (real) Item[3]: 26.273100 (real) Item[4]: 40.769300 (real) Item[5]: 27.212000 (real) Item[6]: 40.192600 (real) Item[7]: 27.210200 (real) Item[3] : GCS Type = dictionary Item[0] : WKT = PROJCS["WGS_1984_UTM_Zone_35N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",50.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",27.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]] (string) Item[1] : Type = PROJCS (name) Item[4] : Type = Measure (name) Item[5] : Bounds Type = array Item[0]: 0 (int) Item[1]: 1 (int) Item[2]: 0 (int) Item[3]: 0 (int) Item[4]: 1 (int) Item[5]: 0 (int) Item[6]: 1 (int) Item[7]: 1 (int) Item[2] : Type = Viewport (name) Item[3] : BBox Type = array Item[0]: 57 (int) Item[1]: 2465 (int) Item[2]: 3052 (int) Item[3]: 46 (int) The resulting GeoTiff has a scaling problem, I suspect related to the transformation matrix. Any thoughts on that? Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] building with PDF driver on Windows + ?
To all, Here's a list of steps for building GDAL with KDE-win's poppler.lib on Windows (win32): 1. Install compiler ( MS VC 9 Express Ed.), or better yet use the full MS VC 9 if you have it. http://www.microsoft.com/express/Downloads/#2008-Visual-CPP I don't think the SDK is required but here's a link: http://blogs.msdn.com/b/windowssdk/archive/2008/02/07/windows-sdk-rtms.aspx 2. Install KDE-win32 stable 4.4.4 (http://windows.kde.org/). Note: In the installer unselect "skip basic settings pages" and select "Package Manager" Select Bin,Devel (and Src?) for: poppler-vc90 0.12.4 freetype-vc902.4.2-3 lcms1.18 3. Get latest GDAL source: svn checkout https://svn.osgeo.org/gdal/trunk/gdal gdal 4. Edit GDAL's nmake.opt to enable poppler (it's commented-outby default), and edit/hack the makefile.vc to disable building of VB6 support (ATL is not included in VC Express Ed. or the SDK) 5. Build GDAL with VC's command line: nmake -f makefile.vc Many thanks to Joaquim and Even for pointing the way... Best Regards, Brent Fraser On 10/17/2010 12:39 PM, Joaquim Luis wrote: On 17-10-2010 19:24, Brent Fraser wrote: Joaquim, Many thanks for the info. I may go the kde-win32 route for now to skip the building of poppler.lib, but I expect that eventually I may need to build it from source, especially for mapserver. Right, and From memory, as I don't find the CMake cash with exact details. obviously I meant "cache" not "cash" ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] building with PDF driver on Windows + ?
Joaquim, Many thanks for the info. I may go the kde-win32 route for now to skip the building of poppler.lib, but I expect that eventually I may need to build it from source, especially for mapserver. Thanks again, Brent Fraser On 10/17/2010 12:11 PM, Joaquim Luis wrote: Brent, From memory, as I don't find the CMake cash with exact details. I first built freetype, which I think is mandatory, and than because I already had them built but don't know if they are a must-be, I added also zlib, jpeg to the CMake list. My first error was to not select 'splash' because I didn't know what it was but it's part of poppler and needs to be selected. I did not include any Qt stuff and probably some other. Once it is built (I used only VC2010) there is the mess of the includes. You will need to include the root directory where the whole poppler tree lives, plus that ../poppler plus still one other dir that holds only one .h file. I moved that .h to one of the other places and now use this (my names) in nmake.opt POPPLER_CFLAGS = -IC:/programs/compa_libs -IC:/programs/compa_libs/poppler Good luck Joaquim I'm about to embark on compiling Poppler on Windows to get Geospatial PDF support in GDAL. Any recommendations with respect to compiler version, dependencies, CMake options, etc? Thanks! Brent Fraser On 10/14/2010 3:23 PM, Joaquim Luis wrote: Did I mention before that the propeller (sorry, poppler) doesn't have any building instructions for Windows? Well, the CMakeLists.txt is incomplete and does not add the contents of the "splash" directory to the project. After adding all *.cc from 'splash' to project, GDAL is happy with the poppler.lib and now gdalinfo says ... KMLSUPEROVERLAY (rw): Kml Super Overlay XYZ (rwv): ASCII Gridded XYZ HF2 (rwv): HF2/HFZ heightfield raster PDF (rov): Geospatial PDF but I'm now confused with the presence of the KMLSUPEROVERLAY because my nmake.opt has # Uncomment out the following to enable KML Super-Overlay support. #KMLSUPEROVERLAY_SUPPORTED = YES #MINIZIP_INCLUDE = -I$(OSSIM_HOME)\minizip\src #MINIZIP_LIBRARY = $(OSSIM_HOME)\minizip\release\minizip.lib Not complaining. Just reporting. Joaquim On 14-10-2010 18:32, Even Rouault wrote: Le jeudi 14 octobre 2010 17:01:06, Joaquim Luis a écrit : Hi, I tried to build GDAL on Win with PDF support and CV2010 Well that's an adventure. I trust you and didn't even try this way. Instead I just downloaded the kde-win32 installer, used the "packager mode" (or whatever it is called. I'm just quoting from my memory of doing this a few weeks ago), and selected the poppler, freetype and lcms packages and their developement packages. The only requirement for GDAL is poppler, freetype and lcms appears to be poppler requirements on this KDE build. But those lib only work with MSVC 2008 if I remember. Even, Hmm, on a further check those symbols are from poppler so it's not a lcms fault. I get another error building poppler referring iconv.h that I ignore either but selected as a no dependency in CMake, but this affect the creation of a poppler-cpp.lib, not poppler.lib so I'm not sure it relates to the GDAL error. I would like to build these dependencies myself because: 1. I don't want to use anything that dares to create manifests dependencies 2. I want to be able to build 64 bits versions too. Anyway, I found these warnings too that are unrelated to this PDF driver issue json_util.c json_util.c(62) : warning C4013: 'open' undefined; assuming extern returning int json_util.c(71) : warning C4013: 'read' undefined; assuming extern returning int json_util.c(74) : warning C4013: 'close' undefined; assuming extern returning int json_util.c(109) : warning C4013: 'write' undefined; assuming extern returning int ogrsf_frmts.lib(resolvexlinks.obj) : warning LNK4221: This object file does not define any previously undefined public symbols, so it will not be used by any link operation that consumes this library and many "blabla ... possible loss of data" warnings. ___ 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] building with PDF driver on Windows + ?
I'm about to embark on compiling Poppler on Windows to get Geospatial PDF support in GDAL. Any recommendations with respect to compiler version, dependencies, CMake options, etc? Thanks! Brent Fraser On 10/14/2010 3:23 PM, Joaquim Luis wrote: Did I mention before that the propeller (sorry, poppler) doesn't have any building instructions for Windows? Well, the CMakeLists.txt is incomplete and does not add the contents of the "splash" directory to the project. After adding all *.cc from 'splash' to project, GDAL is happy with the poppler.lib and now gdalinfo says ... KMLSUPEROVERLAY (rw): Kml Super Overlay XYZ (rwv): ASCII Gridded XYZ HF2 (rwv): HF2/HFZ heightfield raster PDF (rov): Geospatial PDF but I'm now confused with the presence of the KMLSUPEROVERLAY because my nmake.opt has # Uncomment out the following to enable KML Super-Overlay support. #KMLSUPEROVERLAY_SUPPORTED = YES #MINIZIP_INCLUDE = -I$(OSSIM_HOME)\minizip\src #MINIZIP_LIBRARY = $(OSSIM_HOME)\minizip\release\minizip.lib Not complaining. Just reporting. Joaquim On 14-10-2010 18:32, Even Rouault wrote: Le jeudi 14 octobre 2010 17:01:06, Joaquim Luis a écrit : Hi, I tried to build GDAL on Win with PDF support and CV2010 Well that's an adventure. I trust you and didn't even try this way. Instead I just downloaded the kde-win32 installer, used the "packager mode" (or whatever it is called. I'm just quoting from my memory of doing this a few weeks ago), and selected the poppler, freetype and lcms packages and their developement packages. The only requirement for GDAL is poppler, freetype and lcms appears to be poppler requirements on this KDE build. But those lib only work with MSVC 2008 if I remember. Even, Hmm, on a further check those symbols are from poppler so it's not a lcms fault. I get another error building poppler referring iconv.h that I ignore either but selected as a no dependency in CMake, but this affect the creation of a poppler-cpp.lib, not poppler.lib so I'm not sure it relates to the GDAL error. I would like to build these dependencies myself because: 1. I don't want to use anything that dares to create manifests dependencies 2. I want to be able to build 64 bits versions too. Anyway, I found these warnings too that are unrelated to this PDF driver issue json_util.c json_util.c(62) : warning C4013: 'open' undefined; assuming extern returning int json_util.c(71) : warning C4013: 'read' undefined; assuming extern returning int json_util.c(74) : warning C4013: 'close' undefined; assuming extern returning int json_util.c(109) : warning C4013: 'write' undefined; assuming extern returning int ogrsf_frmts.lib(resolvexlinks.obj) : warning LNK4221: This object file does not define any previously undefined public symbols, so it will not be used by any link operation that consumes this library and many "blabla ... possible loss of data" warnings. ___ 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] gdaldem color-relief color_text_files
Mark, Here's a color-by-elevation file from VTP (http://vterrain.org/). You'll have to edit it to get it into GDAL's format: Best Regards, Brent Fraser colormap1 blend: 1 relative: 0 size 35 elev -6800.00 color 255 255 255 elev -6400.00 color 20 20 30 elev -6000.00 color 60 60 70 elev -5600.00 color 120 120 130 elev -5200.00 color 180 185 190 elev -4800.00 color 160 80 0 elev -4400.00 color 128 128 0 elev -4000.00 color 160 0 160 elev -3600.00 color 144 64 144 elev -3200.00 color 128 128 128 elev -2800.00 color 64 128 60 elev -2400.00 color 0 128 0 elev -2000.00 color 0 128 128 elev -1600.00 color 0 0 160 elev -1200.00 color 43 90 142 elev -800.00 color 81 121 172 elev -400.00 color 108 156 195 elev -1.00 color 182 228 255 elev 0.00 color 0 0 238 elev 0.10 color 40 224 40 elev 450.00 color 0 128 0 elev 900.00 color 100 144 76 elev 1350.00 color 204 170 136 elev 1800.00 color 136 100 70 elev 2250.00 color 128 128 128 elev 2700.00 color 180 128 64 elev 3150.00 color 255 144 32 elev 3600.00 color 200 110 80 elev 4050.00 color 160 80 160 elev 4500.00 color 144 40 128 elev 4950.00 color 128 128 128 elev 5400.00 color 255 255 255 elev 5850.00 color 255 255 128 elev 6300.00 color 255 128 0 elev 6750.00 color 0 128 0 On 10/14/2010 6:16 PM, Mark Millman wrote: I am in need of some color_text files to use with the gdaldem color-relief utility to create color relief versions of Slope, Aspect, Roughness, TRI, & TPI maps. I seems to me that these ought to be pretty standard and as my skill set is on the developer side rather than the cartographer side of things I am hoping that some generous cartographer out there has already created some eye candy quality color_text_files for this purpose. Thanks in advance. Mark ___ 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] Working with NED elevation data - part 2
Steve, Excellent info! I had been contemplating using GDAL's VRT format to create an elevation web service so I could display X,Y and Z on mouse moves in web mapping, but was unsure about the performance. Looks like fcgi could work. FYI, I think I'll use Google's XML format (http://code.google.com/apis/maps/documentation/elevation/), but I may extend to have some kind of quality/accuracy and vertical datum values too. Here's a stock Google example: OK 39.7391536 -104.9847034 1608.8402100 36.460 -116.870 -50.7890358 Or perhaps the less verbose USGS format (http://gisdata.usgs.gov/XMLWebServices2/Elevation_Service_Methods.php): [Source Layer ID] [Decimal Elevation Value] [FEET or METERS] Thanks again... Best Regards, Brent Fraser Stephen Woodbridge wrote: Even, It seems the adding: #include resolves the segv issue. I would recommend changing the example in: http://www.gdal.org/gdal_tutorial.html To include that in the section "Reading Raster Data" because the example C code uses CPLMalloc(). It would also be very helpful to add: ... CPLFree(pafScanline); to the example code. And a link to http://gdal.org/cpl__conv_8h.html which is not otherwise referenced. Anyway, code is working great now. I am generating an elevation profile from a number of points along a path. 1. tried to get the points via USGS NED webservice: 132 point took 4-6 min and 132 points is a short path. So I downloaded the USGS NED data and created a .vrt file for the layer: 2. C code to open vrt file and handle requests via local CGI took 4 sec 3. C code to open vrt file and handle requests via local FCGI took < 1 sec And each profile time above has some overhead because it has to compute a route in pgrouting to get the points needed to generate the profile. http://imaptools.com:8080/routeloops/route_loops6.php?x=-71.387764&y=42.620002&length=4828.032&heading=-1&shape=-1&rand=0.06593320066404917&allow=0%2C1&zoom=11&lat=42.62&lon=-71.38776&layers=B00TT&ll=-71.387764%2042.620002&len=3&rs=0.4646901273751023&al=0,1&u=E&cw=1&f=profile Thanks for the help and support everyone. -Steve On 8/22/2010 12:43 PM, Even Rouault wrote: Stephen, at first sight, I don't see anything wrong with your code. You need to include "cpl_conv.h" for CPLMalloc(). And memory allocated with CPLMalloc() should be freed with CPLFree(). The doc is there : http://gdal.org/cpl__conv_8h.html (http://gdal.org/ -> API reference documentation -> Files -> cpl_conv.h ) But that doesn't likely account for the segmentation fault you have. This code should work just fine with a GTiff or VRT dataset. Do you use a recent version of GDAL ? If not, retry ;-) If still crashing, could you attach a debugger and show the stack trace at the place where it crashes and/or run valgrind on it ? Best regards, Even Le dimanche 22 août 2010 05:23:45, Stephen Woodbridge a écrit : OK, well may to answer my own question by reading the GDAL API Tutorial and taking a guess at some parts I came up with the following. This is basically just a copy of the tutorial and at the end I added some code that I'm hoping reads the correct pixel for my lat,lon and prints out the elevation. The tutorial code for "Fetching a Raster Band" works on a GTiff file but segv on a VRT file. Not sure why this happened, although there was a compiler warning about: testit.c: In function 'main': testit.c:80: warning: cast to pointer from integer of different size which would be: float *pafScanline; pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize); I could not find docs on CPLMalloc, or any info on if I need to free that and how! I would appreciate if someone could code review the code at the end and comment on whether or not this is correct. Thanks, -Steve #include int main() { GDALDatasetH hDataset; char *pszFilename = "/u/srcdata/ned2.vrt"; //char *pszFilename = "/u/srcdata/NED2/-90_28_-84_30.tif"; GDALDriverH hDriver; doubleadfGeoTransform[6]; GDALRasterBandH hBand; int nBlockXSize, nBlockYSize; int bGotMin, bGotMax; double adfMinMax[2]; float *pafScanline; int nXSize; int nXOff, nYOff; float elevation; float lon = -74.0; float lat = 42.0; GDALAllRegister(); hDataset = GDALOpen( pszFilename, GA_ReadOnly ); if( hDataset == NULL ) { printf("Failed to open %s!\n", pszFilename); exit(1); } hDriver = GDALGetDatasetDriver( hDataset ); printf( "Driver: %s/%s\n", GDALGetDriverShortName( hDriver ), GDALGetDriverLongName( hDriver ) ); printf( "
Re: [gdal-dev] OGR OVF has no identified FID column
David, GDAL: POINT (0.0 2.0021871) SQL Server: POINT (150 150) That IS strange! There could be a problem with GDAL's or SQL Server's encoding/decoding of WKB, or it could be related to data type (geometry? geography?) or the lack of a Spatial Reference System? Best Regards, Brent Fraser David Lowther wrote: D:\Utility\ms4w\tools\gdal-ogr>ogrinfo test.ovf PointData --debug on ERROR 4: Update access not supported for VRT datasources. OGR_ODBC: EstablishSession(DSN:"CLO", userid:"un", password:"pw") ODBC: SQLConnect(CLO) OGR_ODBC: Table PointData has no identified FID column. OGR: OGROpen(ODBC:un/p...@clo,PointData/00FB9BB0) succeeded as ODBC. OGR: OGROpen(test.ovf/010870B8) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 Layer name: PointData Geometry: Unknown (any) OGR_ODBC: ExecuteSQL(SELECT * FROM PointData) Feature Count: 1 OGR_ODBC: ExecuteSQL(SELECT * FROM PointData) Extent: (0.00, 2.00) - (0.00, 2.00) Layer SRS WKT: (unknown) PKey: Integer (10.0) Feature: Binary (0.0) OGR_ODBC: ExecuteSQL(SELECT * FROM PointData) OGRFeature(PointData):1 PKey (Integer) = 1 Feature (Binary) = 010C00C0624000C06240 POINT (0.0 2.0021871) VRT: 2 features read on layer 'PointData'. OGR_ODBC: 3 features read on layer 'PointData'. ODBC: SQLDisconnect() And, why do you think it is saying 0.0 2.00... for the point when the SQL query says? 1 0x010C00C0624000C06240 POINT (150 150) David Lowther Coordinate Solutions, Inc. -Original Message- From: Brent Fraser [mailto:bfra...@geoanalytic.com] Sent: Thursday, August 19, 2010 5:21 PM To: David Lowther Cc: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] OGR OVF has no identified FID column I think it's faster since it doesn't have to enumerate all the tables,views, etc. so how do you get for attributes names when you do a ogrinfo test.ovf PointData --debug on Brent David Lowther wrote: Still doesn't see the FID but is way faster! ODBC:un/p...@clo,PointData PKey ogrinfo test.ovf --debug on ERROR 4: Update access not supported for VRT datasources. OGR_ODBC: EstablishSession(DSN:"CLO", userid:"un", password:"pw") ODBC: SQLConnect(CLO) OGR_ODBC: Table PointData has no identified FID column. OGR: OGROpen(ODBC:un/p...@clo,PointData/02799B88) succeeded as ODBC. OGR: OGROpen(test.ovf/010270B8) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 1: PointData ODBC: SQLDisconnect() Does the datatype of the FID column matter? I have it as an int, does it need to be guid? David Lowther Coordinate Solutions, Inc. -Original Message- From: Brent Fraser [mailto:bfra...@geoanalytic.com] Sent: Thursday, August 19, 2010 3:49 PM To: David Lowther Cc: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] OGR OVF has no identified FID column Yeah I saw that after I sent the message. Take out the , it's really only useful if want to do "special" things, and it can kill spatial indexing on ODBC datasources. Does this work any better: ODBC:un/p...@clo,PointData PKey Brent David Lowther wrote: Brent, Thanks for the reply. On initial read I thought you had found my issue - that I was trying so many things I finally messed it all up. But on closer inspection I am aliasing the PKey column to FID in the SrcSQL (select PKey as FID, Feature from PointData). When FID ogrinfo returns: ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. OGR: OGROpen(test.ovf/00ED7090) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 1: UserPoint OGR: ReleaseDataSource(ODBC:un/p...@clo/02719D10) dereferenced and now destroying. ODBC: SQLDisconnect() When PKEY ogrinfo returns: ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. ERROR 1: Unable to identify FID field 'PKey'. OGR: ReleaseDataSource(ODBC:un/p...@clo/02C39D18) dereferenced and now destroying. ODBC: SQLDisconnect() FAILURE: Unable to open datasource `test.ovf' with the following drivers. -> ESRI Shapefile David Lowther Coordinate Solutions, Inc. -Original Message- From: Brent Fraser [mailto:bfra...@geoanalytic.com] Sent: Thursday, August 19, 2010 3:15 PM To: David Lowther Cc: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] OGR OVF has no identified FID column How about using: PKey but that doesn
Re: [gdal-dev] OGR OVF has no identified FID column
I think it's faster since it doesn't have to enumerate all the tables,views, etc. so how do you get for attributes names when you do a ogrinfo test.ovf PointData --debug on Brent David Lowther wrote: Still doesn't see the FID but is way faster! ODBC:un/p...@clo,PointData PKey ogrinfo test.ovf --debug on ERROR 4: Update access not supported for VRT datasources. OGR_ODBC: EstablishSession(DSN:"CLO", userid:"un", password:"pw") ODBC: SQLConnect(CLO) OGR_ODBC: Table PointData has no identified FID column. OGR: OGROpen(ODBC:un/p...@clo,PointData/02799B88) succeeded as ODBC. OGR: OGROpen(test.ovf/010270B8) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 1: PointData ODBC: SQLDisconnect() Does the datatype of the FID column matter? I have it as an int, does it need to be guid? David Lowther Coordinate Solutions, Inc. -Original Message- From: Brent Fraser [mailto:bfra...@geoanalytic.com] Sent: Thursday, August 19, 2010 3:49 PM To: David Lowther Cc: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] OGR OVF has no identified FID column Yeah I saw that after I sent the message. Take out the , it's really only useful if want to do "special" things, and it can kill spatial indexing on ODBC datasources. Does this work any better: ODBC:un/p...@clo,PointData PKey Brent David Lowther wrote: Brent, Thanks for the reply. On initial read I thought you had found my issue - that I was trying so many things I finally messed it all up. But on closer inspection I am aliasing the PKey column to FID in the SrcSQL (select PKey as FID, Feature from PointData). When FID ogrinfo returns: ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. OGR: OGROpen(test.ovf/00ED7090) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 1: UserPoint OGR: ReleaseDataSource(ODBC:un/p...@clo/02719D10) dereferenced and now destroying. ODBC: SQLDisconnect() When PKEY ogrinfo returns: ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. ERROR 1: Unable to identify FID field 'PKey'. OGR: ReleaseDataSource(ODBC:un/p...@clo/02C39D18) dereferenced and now destroying. ODBC: SQLDisconnect() FAILURE: Unable to open datasource `test.ovf' with the following drivers. -> ESRI Shapefile David Lowther Coordinate Solutions, Inc. -Original Message- From: Brent Fraser [mailto:bfra...@geoanalytic.com] Sent: Thursday, August 19, 2010 3:15 PM To: David Lowther Cc: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] OGR OVF has no identified FID column How about using: PKey but that doesn't explain why it takes a long time... Best Regards, Brent Fraser David Lowther wrote: List, I have an OVF layer defined as follows: ODBC:un/p...@clo select PKey as FID, Feature from PointData FID The ODBC connection is to a SQL Express 2008 database. The PointData table is defined as follows: CREATE TABLE [dbo].[PointData]( [PKey] [int] IDENTITY(1,1) NOT NULL, [Feature] [geometry] NULL, CONSTRAINT [PK_PointData] PRIMARY KEY CLUSTERED ( [PKey] ASC )WITH (PAD_INDEX = OFF, STATISTICS_NORECOMPUTE = OFF, IGNORE_DUP_KEY = OFF, ALLOW_ROW_LOCKS = ON, ALLOW_PAGE_LOCKS = ON) ON [PRIMARY] ) ON [PRIMARY] This query select *, Feature.STAsText() from pointdata yields 1 0x010C00C0624000C06240 POINT (150 150) "ogrinfo test.ovf --debug on" runs for a very long time returning many lines of this sort of stuff: ERROR 4: Update access not supported for VRT datasources. OGR_ODBC: EstablishSession(DSN:"CLO", userid:"un", password:"pw") ODBC: SQLConnect(CLO) ODBC: CatalogNameL: (null) Schema name: (null) OGR_ODBC: Table PointData has no identified FID column. OGR_ODBC: Table UserMapData has no identified FID column. OGR_ODBC: Table CHECK_CONSTRAINTS has no identified FID column. before finally saying: OGR: OGROpen(ODBC:un/p...@clo/02C79D10) succeeded as ODBC. ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. OGR: OGROpen(test.ovf/00E17090) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 1: UserPoint OGR: ReleaseDataSource(ODBC:un/p...@clo/02C79D10) dereferenced and now destroying. ODBC: SQLDisconnect() GDAL version is 1.6.0. I have tried using no FID tag and using SrcLayer instead of SrcSQL. Anybody know what I could do
Re: [gdal-dev] OGR OVF has no identified FID column
Yeah I saw that after I sent the message. Take out the , it's really only useful if want to do "special" things, and it can kill spatial indexing on ODBC datasources. Does this work any better: ODBC:un/p...@clo,PointData PKey Brent David Lowther wrote: Brent, Thanks for the reply. On initial read I thought you had found my issue - that I was trying so many things I finally messed it all up. But on closer inspection I am aliasing the PKey column to FID in the SrcSQL (select PKey as FID, Feature from PointData). When FID ogrinfo returns: ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. OGR: OGROpen(test.ovf/00ED7090) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 1: UserPoint OGR: ReleaseDataSource(ODBC:un/p...@clo/02719D10) dereferenced and now destroying. ODBC: SQLDisconnect() When PKEY ogrinfo returns: ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. ERROR 1: Unable to identify FID field 'PKey'. OGR: ReleaseDataSource(ODBC:un/p...@clo/02C39D18) dereferenced and now destroying. ODBC: SQLDisconnect() FAILURE: Unable to open datasource `test.ovf' with the following drivers. -> ESRI Shapefile David Lowther Coordinate Solutions, Inc. -----Original Message- From: Brent Fraser [mailto:bfra...@geoanalytic.com] Sent: Thursday, August 19, 2010 3:15 PM To: David Lowther Cc: gdal-dev@lists.osgeo.org Subject: Re: [gdal-dev] OGR OVF has no identified FID column How about using: PKey but that doesn't explain why it takes a long time... Best Regards, Brent Fraser David Lowther wrote: List, I have an OVF layer defined as follows: ODBC:un/p...@clo select PKey as FID, Feature from PointData FID The ODBC connection is to a SQL Express 2008 database. The PointData table is defined as follows: CREATE TABLE [dbo].[PointData]( [PKey] [int] IDENTITY(1,1) NOT NULL, [Feature] [geometry] NULL, CONSTRAINT [PK_PointData] PRIMARY KEY CLUSTERED ( [PKey] ASC )WITH (PAD_INDEX = OFF, STATISTICS_NORECOMPUTE = OFF, IGNORE_DUP_KEY = OFF, ALLOW_ROW_LOCKS = ON, ALLOW_PAGE_LOCKS = ON) ON [PRIMARY] ) ON [PRIMARY] This query select *, Feature.STAsText() from pointdata yields 1 0x010C00C0624000C06240 POINT (150 150) "ogrinfo test.ovf --debug on" runs for a very long time returning many lines of this sort of stuff: ERROR 4: Update access not supported for VRT datasources. OGR_ODBC: EstablishSession(DSN:"CLO", userid:"un", password:"pw") ODBC: SQLConnect(CLO) ODBC: CatalogNameL: (null) Schema name: (null) OGR_ODBC: Table PointData has no identified FID column. OGR_ODBC: Table UserMapData has no identified FID column. OGR_ODBC: Table CHECK_CONSTRAINTS has no identified FID column. before finally saying: OGR: OGROpen(ODBC:un/p...@clo/02C79D10) succeeded as ODBC. ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called. OGR_ODBC: Table SELECT has no identified FID column. OGR: OGROpen(test.ovf/00E17090) succeeded as VRT. Had to open data source read-only. INFO: Open of `test.ovf' using driver `VRT' successful. OGR: GetLayerCount() = 1 1: UserPoint OGR: ReleaseDataSource(ODBC:un/p...@clo/02C79D10) dereferenced and now destroying. ODBC: SQLDisconnect() GDAL version is 1.6.0. I have tried using no FID tag and using SrcLayer instead of SrcSQL. Anybody know what I could do to make it see the FID? I'm guessing it would be a lot faster if it didn't check every table in the database for an FID... Thanks in advance for your assistance, Dave Lowther ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Raster Strategies Examples
The steps you've listed below are correct for maximizing performance when viewing the raster at (or around) the native resolution of the data. However, as you've found out, the performance can still be poor when trying to view the entire dataset. The basic rule is that your application (mapserver?) should not have to open more than 4 to 8 raster files to render the view. With that in mind (and since disk space is not a problem), create "layers" of the raster dataset at different scales, say a factor of 2, until you have only one image. For example, if the resolution of your original data is 1 meter per pixel, your layers would look like the following. Note the first eight layers are handled by the internal overviews you created with gdaladdo. Pixel Map Scale: Size: 1:n - -- 1 m 4,000 2 8,000 4 16,000 8 32k 16 64 32 128 64 250 128 500 Additional External Layers 256 1m 512 2m 10244m 20488m 409616m To create the additional external layers, use GDAL's gdalbuiltvrt to create a VRT file of the original dataset. You can then use gdal_retile.py or Maptiler to build the layers. Best Regards, Brent Fraser Edi KARADUMI wrote: I have read many strategies for raster performance, but i still have problems with my case. The posts that i have read explain the strategie, but are not very detailed. Im new to mapserver so i have problems implementing them. My case is: -about 6000 tiles that form the map -170mb each tile -aproximately 1.2T of iamges -i have 10T HD and i think disk space is not a problem for me -tile format is TIF -tile size 9375x6250 pixels The strategy i have implemented and the problems i have -First i divided the tiles in 60 folders, to increase the performance in disk seek/read. -Than i transformed source files into internally tiles with the command gdal_translate -of GTiff -co "TILED=YES" original.tif tiled.tif -i used gdaladdo to add internal pyramids gdaladdo -r average 2 4 8 16 32 64 128 -Created a tileindex using gdaltindex -write_absolute_path MapAll.shp //server/Maps/Subfolder1/*.tif -Created a spatial index .qix file shptree MapAll.shp -than added the layer to the mapfile without the .shp extension so the application can use the .qix as you may know, i have very slow performance when i zoom out and im stucked here. As i have read i should make a copy of the tiles with reduced resolution. Merge the tiles together and use min/max scale to show different layers in different scales. the min/max scales i zoom in/out are 100/120. Now my questions are - How can i calculate the scale where i should create another layer of the tiles, or i shoud see it with some tests? - how much should be the resolution of the new layer? - is there any tools or program to merge the tiles? merging 6000 tiles with the gdalwarp by writing the command by myself is frustrating - how many tiles should i merge together to create the new layer? (how many tiles should have the new layer) i know that in each zoomscale its better to appear only one tile but i dont know how to calculate it - the tiles that i should merge are the originals or those with internal tiling and overviews? ___ 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] ERROR 1: Too many points (441 out of 441) failed to transform, , unable to compute output bounds.
Stephan, Even though the tiles are South of the equator, they are in the UTM North coordinate system. Try -a_srs EPSG:32601 instead of 32701, for Zone 1. Brent Stephen Woodbridge wrote: Even or anyone, I have made some progress, I think, in that all of the northern hemisphere files have been translated without error, but until I can get it working in mapserver to assess the results, I wont be sure. But all of the southern hemisphere files now generate this error regardless of using --config CENTER_LONG 180 or not. So what is the "standard" if there is one way of assessing how to deal with files like this. I have 882 .sid files 286 in the south and 596 in the north. I can script the conversion, but I'm not sure what I need to look for and what I need to do for any random file, like setting --config CENTER_LONG 180 or doing something else. It seems like having a tool like this would make sense and be useful to others. So here is the info on the first of the southern files: + rm -f /geotiff/GeoCover2000TIF//s-01-00.tif + gdal_translate -a_srs EPSG:32701 -of VRT /geotiff/GeoCover2000/s-01-00.sid -b 1 -b 2 -b 3 /geotiff/GeoCover2000TIF/tmp.vrt Input file size is 26981, 38945 + gdalwarp -srcnodata 0 -dstnodata 0 -s_srs EPSG:32701 -t_srs EPSG:4326 -rb -wm 250 --config GDAL_CACHEMAX 250 --config GDAL_ONE_BIG_READ ON -co TILED=YES /geotiff/GeoCover2000TIF/tmp.vrt /geotiff/GeoCover2000TIF//s-01-00.tif ERROR 1: Too many points (440 out of 441) failed to transform, unable to compute output bounds. gdalinfo /geotiff/GeoCover2000TIF/tmp.vrt Driver: VRT/Virtual Raster Files: /geotiff/GeoCover2000TIF/tmp.vrt /geotiff/GeoCover2000/s-01-00.sid Size is 26981, 38945 Coordinate System is: PROJCS["WGS 84 / UTM zone 1S", 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"]], UNIT["metre",1, AUTHORITY["EPSG","9001"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-177], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",50], PARAMETER["false_northing",1000], AUTHORITY["EPSG","32701"], AXIS["Easting",EAST], AXIS["Northing",NORTH]] Origin = (499498.125,1011.750) Pixel Size = (14.250,-14.250) Metadata: IMAGE__INPUT_NAME=/dbuy/l7/outputs/mosaics742/S-01-0/locals/S-01-00_lr_742_2000.tif, /dbuy/l7/outputs/mosaics742/S-01-0/locals/S-01-00_ur_742_2000.tif IMAGE__INPUT_FILE_SIZE=3164295044.00 IMAGE__COMPRESSION_VERSION=1,6,3 IMAGE__TARGET_COMPRESSION_RATIO=29.98 IMAGE__COMPRESSION_NLEV=8 IMAGE__COMPRESSION_WEIGHT=2.00 IMAGE__COMPRESSION_GAMMA=1.00 IMAGE__COMPRESSION_BLOCK_SIZE=4096 IMAGE__CREATION_DATE=Mon Sep 22 13:29:44 2003 VERSION=MG2 Corner Coordinates: Upper Left ( 499498.125,1011.750) (177d 0'0.00"W, 90d 0'0.00"S) Lower Left ( 499498.125, -553954.500) (177d 0'0.00"W, 90d 0'0.00"S) Upper Right ( 883977.375,1011.750) (177d 0'0.00"W, 90d 0'0.00"S) Lower Right ( 883977.375, -553954.500) (177d 0'0.00"W, 90d 0'0.00"S) Center ( 691737.750, -276471.375) (177d 0'0.00"W, 90d 0'0.00"S) Band 1 Block=128x128 Type=Byte, ColorInterp=Red Band 2 Block=128x128 Type=Byte, ColorInterp=Green Band 3 Block=128x128 Type=Byte, ColorInterp=Blue Thanks, -Steve Even Rouault wrote: So does --config CENTER_LONG 180 impact processing of files that are in other UTM Zones? ie: can I just always use this or will if mess up other files in some way? No, I'd advise against always using --config CENTER_LONG 180 if you want to be able to deal with what lies at the west of the Greenwich meridian... (That would probably work but the longitudes would not fit in the [-180,180] interval as expected) Here's why (ogr/ogrct.cpp) : if( x[i] < dfTargetWrapLong - 180.0 ) x[i] += 360.0; else if( x[i] > dfTargetWrapLong + 180 ) x[i] -= 360.0; } So, if the coordinates are too far (more than 180 deg) from the center longitude, 360 is added or substracted. That helps when you're warping from UTM 1 to EPSG:4326 as part of the target coordinates are > 179 and part < -179. The +360 will correct the ones < -179 so that they end up being at the east of the ones > 179. I imagine that some heuristics could be added to gdalwarp to detect your situation and automatically add the --config CENTER_LONG 180, but there's always the risk to break other working use cases, and anyway people should be aware that there's probably some post-processing needed to deal with the fact that the result doesn't fit
Re: [gdal-dev] OGR and SpatiaLite
Thanks! I keep forgetting about those builds. But I do think the SQLite driver should report if it supports SpatiaLite. For example, in ogrsqlitedriver.cpp have something like: const char *OGRSQLiteDriver::GetName() { #ifdef HAVE_SPATIALITE return "SQLite+SpatiaLite functions"; #else return "SQLite"; #endif } Although looking at the driver code, the support really isn't as simple as that. As the OGR driver doc mentions the OGR implementation of SQLite will read and write SpatiaLite databases, it just will not do spatial indexes or SpatiaLite spatial functions (as these are in the code library, not in the database as Stored Procedures as with other spatial database types like PostGIS). But it WILL read and write SpatiaLite's flavour of geometry, and create SpatiaLite-compatible databases (it even creates the geometry_columns and spatial_ref_sys tables!). I think these points need to be stated explicitly in the driver doc. Best Regards, Brent Fraser Tamas Szekeres wrote: Check out the recent builds from here: http://vbkto.dyndns.org/sdk/ Those have been compiled with spatialite support. Best regards, Tamas 2010/3/11 Brent Fraser <mailto:bfra...@geoanalytic.com>> Tamas, I was afraid of that. I'm having trouble getting MS4W 3 Beta 10's OGR to read my SpatiaLite database. I suspect it's because SpatiaLite was not compiled in (I'm a SQLite newbie so I'm not sure); but there's no way to tell. Thanks! Brent > Brent, > > It should be "SQLite" in both cases. > > Best regards, > > Tamas > > > 2010/3/11 Brent Fraser mailto:bfra...@geoanalytic.com>> > > > What would "ogrinfo --formats" show if GDAL was compiled with > > HAVE_SPATIALITE set? > > What would it show if only SQLite (and not SpatiaLite) was compiled in? > > > > Thanks! > > Brent Fraser > > > > > > ___ > > gdal-dev mailing list > > gdal-dev@lists.osgeo.org <mailto: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] OGR and SpatiaLite
Execellent! Thanks! Brent Yewondwossen Assefa wrote: Brent, MS4W 3 Beta 10's OGR was not build with that -DHAVE_SPATIALITE. Next version would have it (bug refernece: http://bugzilla.maptools.org/show_bug.cgi?id=2167) regards, Brent Fraser wrote: Tamas, I was afraid of that. I'm having trouble getting MS4W 3 Beta 10's OGR to read my SpatiaLite database. I suspect it's because SpatiaLite was not compiled in (I'm a SQLite newbie so I'm not sure); but there's no way to tell. Thanks! Brent Brent, It should be "SQLite" in both cases. Best regards, Tamas 2010/3/11 Brent Fraser What would "ogrinfo --formats" show if GDAL was compiled with HAVE_SPATIALITE set? What would it show if only SQLite (and not SpatiaLite) was compiled in? Thanks! Brent Fraser ___ 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] OGR and SpatiaLite
Tamas, I was afraid of that. I'm having trouble getting MS4W 3 Beta 10's OGR to read my SpatiaLite database. I suspect it's because SpatiaLite was not compiled in (I'm a SQLite newbie so I'm not sure); but there's no way to tell. Thanks! Brent > Brent, > > It should be "SQLite" in both cases. > > Best regards, > > Tamas > > > 2010/3/11 Brent Fraser > >> What would "ogrinfo --formats" show if GDAL was compiled with >> HAVE_SPATIALITE set? >> What would it show if only SQLite (and not SpatiaLite) was compiled in? >> >> Thanks! >> Brent Fraser >> >> >> ___ >> 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] OGR and SpatiaLite
What would "ogrinfo --formats" show if GDAL was compiled with HAVE_SPATIALITE set? What would it show if only SQLite (and not SpatiaLite) was compiled in? Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] extract vector/raster data from GeoPDF
User error. I had to turn on the Tools -> Analysis -> Geospatial Location Tool to get the scrolling Lat/Lon. Brent Fraser wrote: Interesting. I downloaded it with Firefox, opened it in Acrobat 9, but it doesn't recognize it as a GeoPDF. Hmmm... Joaquim Luis wrote: Brent Fraser wrote: Since the TerraGo way has been superceded by the Adobe/ISO way (thank goodness!), it makes it easier to identify the reference system definition in GDAL, so there may be some hope for GDAL and GeoPDF. Still on this matter, the GeoPDF http://www.terragotech.com/blog/quada.pdf written in the Adobe way is not recognized as a GeoPDF by my Acrobat with the TerraGo Toolbar and it crashes Chrome when trying to download it!!! ___ 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] extract vector/raster data from GeoPDF
Interesting. I downloaded it with Firefox, opened it in Acrobat 9, but it doesn't recognize it as a GeoPDF. Hmmm... Joaquim Luis wrote: Brent Fraser wrote: Since the TerraGo way has been superceded by the Adobe/ISO way (thank goodness!), it makes it easier to identify the reference system definition in GDAL, so there may be some hope for GDAL and GeoPDF. Still on this matter, the GeoPDF http://www.terragotech.com/blog/quada.pdf written in the Adobe way is not recognized as a GeoPDF by my Acrobat with the TerraGo Toolbar and it crashes Chrome when trying to download it!!! ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] extract vector/raster data from GeoPDF
Joaquim, Interesting! As Peter pointed out, the LPTS tag is the transformation matrix to go from projected coords to page coords, and is created by forming the UTM->page coordinates tranformation (from the information in the CTM tag below) and inverting it (there's a brief mention near the end of http://geopdf.blogspot.com/2009/02/georegistration-worked-example.html) Since the TerraGo way has been superceded by the Adobe/ISO way (thank goodness!), it makes it easier to identify the reference system definition in GDAL, so there may be some hope for GDAL and GeoPDF. Now I need to think about my requirements for reading and writing GeoPDFs. Write a Postscript driver for GDAL? Figure out how to write the equivalent info into a PDF using Mapserver/Cairo? Yikes! Thanks, Brent Joaquim Luis wrote: Brent, I was also digging a bit on this matter. It's true that the blogger only opens our appetite as there is no more details on how to do it generally. I mean, for UTM maps it has become fairly obvious how one can do it. The relevant part in the PS file is [ {ThisPage} << /LGIDict << /Description Title /CTM [(35.28267) (0) (0) (35.28267) (205188.64) (3207094.8)] /Projection << /Description (WGS 84 UTM 16N) /ProjectionType (TC) /Datum (WGE) /CentralMeridian (-87.0) /OriginLatitude (0.0) /FalseEasting (50.0) /FalseNorthing (0.0) /ScaleFactor (0.999600) /Type /Projection >> /Type /LGIDict /Version (2.1) > > /PUT pdfmark However, there is no more info on how to program in this "TerraGo way" A next post shads a bit more lite on the matter. http://geopdf.blogspot.com/2009/02/geopdf-and-geops-with-adobe-style.html As you can see in the quada.ps file there is an alternative "Adobe way", which encodes the referencing as % embed georegistation info [ {ThisPage} << /VP [ << /Type /Viewport /BBox[0 0 dh dv] /Name Title /Measure << /Type /Measure /Subtype /GEO /Bounds[0 0 0 1 1 1 1 0] /GPTS[29.0 -90.0 30.0 -90.0 30.0 -89.0 29.0 -89.0] /LPTS[0.024324 0.039461 0.051689 0.978774 0.975675 0.961136 0.957432 0.022122] /GCS << /Type /PROJCS /WKT (PROJCS["WGS_1984_UTM_Zone_16N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",50],PARAMETER["false_northing",0],UNIT["Meter",1]]) >> >> >>] > /PUT pdfmark Now this is more familiar, except for the /LPTS[0.024324 0.039461 0.051689 0.978774 0.975675 0.961136 0.957432 0.022122] part that, contrary to the promise has not been explained where it comes from. I think I could implement the first encoding type in GMT's ps2raster (if I knew all the specs) Second type would be even easier ... if we had a way to translate the GMT projection syntax into WKT, which we don't. Joaquim Joaquim, Doing the format conversion from .ps to pdf is one thing (and there are several ways to do it), but embedding the georeferencing in the PDF to make it a GeoPDF is the interesting bit. I haven't found any open source project capable of do that. Or even capable of writing the PDF elements (e.g. Frames) necessary for writing the GeoPDF information. But I remain hopeful... Brent Joaquim Luis wrote: Brent Fraser wrote: I see the example calculation of the transformation matrix, and a statement "I created a GeoPDF by running the Postscript file through Ghostscript to create a Postscript file that looks like this" But I don't see where the georeferencing is written into PDF to make it a GeoPDF. Is it a Ghostscript command line? Brent, One only has to convert it to pdf using ghostscript. I did it with GMT's ps2raster, like that ps2raster quad.ps -Tf -A Joaquim Luis Brent Klokan Petr Přidal wrote: Hi, There is a great blog post (and the linked "worked example" post with details): http://geopdf.blogspot.com/2009/02/geopdf-map-for-worked-example.html It shows you how to create geopdf via GhostScript - so there is already a practical open-source example how to encode the georeference into the PDF/PS according the OGC standard - for use in Acrobat Reader. To add support for such tag in MapServer, which generates pdf dynamicaly via pdflib, should not be totally problematic. Decoding is not as hard either, there are nice libraries like poppler (http://poppler.freedesktop.org/), which allows you to parse vectors (and convert them to SVG for example) or rasterize the PDF file
Re: [gdal-dev] extract vector/raster data from GeoPDF
Joaquim, Doing the format conversion from .ps to pdf is one thing (and there are several ways to do it), but embedding the georeferencing in the PDF to make it a GeoPDF is the interesting bit. I haven't found any open source project capable of do that. Or even capable of writing the PDF elements (e.g. Frames) necessary for writing the GeoPDF information. But I remain hopeful... Brent Joaquim Luis wrote: Brent Fraser wrote: I see the example calculation of the transformation matrix, and a statement "I created a GeoPDF by running the Postscript file through Ghostscript to create a Postscript file that looks like this" But I don't see where the georeferencing is written into PDF to make it a GeoPDF. Is it a Ghostscript command line? Brent, One only has to convert it to pdf using ghostscript. I did it with GMT's ps2raster, like that ps2raster quad.ps -Tf -A Joaquim Luis Brent Klokan Petr Přidal wrote: Hi, There is a great blog post (and the linked "worked example" post with details): http://geopdf.blogspot.com/2009/02/geopdf-map-for-worked-example.html It shows you how to create geopdf via GhostScript - so there is already a practical open-source example how to encode the georeference into the PDF/PS according the OGC standard - for use in Acrobat Reader. To add support for such tag in MapServer, which generates pdf dynamicaly via pdflib, should not be totally problematic. Decoding is not as hard either, there are nice libraries like poppler (http://poppler.freedesktop.org/), which allows you to parse vectors (and convert them to SVG for example) or rasterize the PDF files (into TIFF,...) via Cairo. The work is in assigning correct geographic coordinates to the coordinate system internally used in PDF files and especially write the bridge to the outside world (with GDAL/OGR). I am afraid that authors of the GeoPDF standard would not like this, as it seems that the idea of GeoPDF is "see it in the Acrobat, print it, but that's all". At least I think so, because they discontinued their Geopdf2geotiff product and all the conversion tools are just one way - into GeoPDF. Please correct me... Anyway, in this moment you can quite easily use utility like "pdfimages" to extract full quality image tiles from any GeoPDF (like those from USGS) and merge it based on their location in PDF into one GDAL file via VRT (gdalbuildvrt) with a bit of hacking. This is what I did for my favorite USGS DRG of Grand Canyon ;-). Look at: http://klokan.mzk.cz/~klokan/geopdf/ - soon I will update the MapTiler.org overlay examples... Unfortunately all PDF parsing libraries I know are GPL, and that means we can't use them for the gdal driver - because of the license issues. But to create a GPL utility for converting GeoPDF to anything what GDAL/OGR supports should be OK. Poppler can be the best base of such GDAL-based utility for reading/rasterizing of the GeoPDF files. Now just find a sponsor and time to make it ;-). Best, Klokan Petr Pridal ___ 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] extract vector/raster data from GeoPDF
I see the example calculation of the transformation matrix, and a statement "I created a GeoPDF by running the Postscript file through Ghostscript to create a Postscript file that looks like this" But I don't see where the georeferencing is written into PDF to make it a GeoPDF. Is it a Ghostscript command line? Brent Klokan Petr Přidal wrote: Hi, There is a great blog post (and the linked "worked example" post with details): http://geopdf.blogspot.com/2009/02/geopdf-map-for-worked-example.html It shows you how to create geopdf via GhostScript - so there is already a practical open-source example how to encode the georeference into the PDF/PS according the OGC standard - for use in Acrobat Reader. To add support for such tag in MapServer, which generates pdf dynamicaly via pdflib, should not be totally problematic. Decoding is not as hard either, there are nice libraries like poppler (http://poppler.freedesktop.org/), which allows you to parse vectors (and convert them to SVG for example) or rasterize the PDF files (into TIFF,...) via Cairo. The work is in assigning correct geographic coordinates to the coordinate system internally used in PDF files and especially write the bridge to the outside world (with GDAL/OGR). I am afraid that authors of the GeoPDF standard would not like this, as it seems that the idea of GeoPDF is "see it in the Acrobat, print it, but that's all". At least I think so, because they discontinued their Geopdf2geotiff product and all the conversion tools are just one way - into GeoPDF. Please correct me... Anyway, in this moment you can quite easily use utility like "pdfimages" to extract full quality image tiles from any GeoPDF (like those from USGS) and merge it based on their location in PDF into one GDAL file via VRT (gdalbuildvrt) with a bit of hacking. This is what I did for my favorite USGS DRG of Grand Canyon ;-). Look at: http://klokan.mzk.cz/~klokan/geopdf/ - soon I will update the MapTiler.org overlay examples... Unfortunately all PDF parsing libraries I know are GPL, and that means we can't use them for the gdal driver - because of the license issues. But to create a GPL utility for converting GeoPDF to anything what GDAL/OGR supports should be OK. Poppler can be the best base of such GDAL-based utility for reading/rasterizing of the GeoPDF files. Now just find a sponsor and time to make it ;-). Best, Klokan Petr Pridal ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] extract vector/raster data from GeoPDF
Keith, No talk that I'm aware of. There was a little chatter on the Geowanking list: http://geowanking.org/pipermail/geowanking_geowanking.org/2009-February/017492.html And occasionally some interest about producing PDF from Mapserver (e.g. http://lists.osgeo.org/pipermail/mapserver-users/2009-July/062035.html). I'm hoping that as GeoPDF becomes popular some funding can be found to build a GDAL/OGR driver to read it. Best Regards, Brent Fraser keith lewis wrote: Hello people, Is there any talk about developing tools in GDAL/OGR to extract raster or vector data from a GeoPDF now the 2.2 version was submitted to the OGC?? The submitted specs are for the georeferencing aspect only. ArcGIS and TerraGo tools can create a GeoPDF. The USGS offers free Geopdf that contain geospatial data. I'm sure someone will eventually release a tool. Keith ___ 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] OGR and database filter on spatial extent
Frank, That did it! I was able to get my query back down to < 2 seconds using: SELECT * FROM alta83_v1 WHERE BOT_LONG > -114 AND BOT_LONG < -113 AND BOT_LAT > 51 AND BOT_LAT < 52 Now I've got to enhance GeoMoose to calculate the view extents in Lat/Long and hope mapserver will like a VRT layer def of something like: SELECT * FROM alta83_v1 WHERE BOT_LONG > %VIEW_MIN_LONG% AND BOT_LONG < %VIEW_MAX_LONG% AND BOT_LAT > %VIEW_MIN_LAT% AND BOT_LAT < %VIEW_MAX_LAT% Many thanks! Brent Frank Warmerdam wrote: Brent Fraser wrote: Frank, > Brent, Any time you provide your own SQL it will be executed unaltered. The spatial and attribute filters are only automatically built into the select statement when the select statement is constructed internally by OGR. OGR makes *no attempt at all* to interprete provided SQL (via ExecuteSQL() or . So it has no idea what tables it relates to or how any restrictions could be added. I guess that's technically true, but ogr (likely the VRT parser) gives an error when I try to construct my own spatial filter: SELECT * FROM alta83_v1 WHERE BOT_LONG > -114 AND BOT_LONG < -113 AND BOT_LAT > 51 AND BOT_LAT < 52 ERROR 1: Line 4: Didn't find expected '=' for value of attribute 'AND'. I expect it doesn't like the ">". Quoting the SQL and bracketing the where clause portion doesn't help. Perhaps there's an escape sequence I could use? Brent, You need to use normal XML escaping. I imagine that is > and < for the comparison operators. And is there point in constructing a view with columns of WKT geometry, XMIN, YMIN, XMAX and YMAX to see if the ODBC driver would take advantage of those? Not if you are going to use SrcSQL - no. If you provide your own SQL, then no automatic spatial filtering will take place in the database - it will be done after fetching all rows. Best regards, ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR and database filter on spatial extent
Frank, > Brent, Any time you provide your own SQL it will be executed unaltered. The spatial and attribute filters are only automatically built into the select statement when the select statement is constructed internally by OGR. OGR makes *no attempt at all* to interprete provided SQL (via ExecuteSQL() or . So it has no idea what tables it relates to or how any restrictions could be added. I guess that's technically true, but ogr (likely the VRT parser) gives an error when I try to construct my own spatial filter: SELECT * FROM alta83_v1 WHERE BOT_LONG > -114 AND BOT_LONG < -113 AND BOT_LAT > 51 AND BOT_LAT < 52 ERROR 1: Line 4: Didn't find expected '=' for value of attribute 'AND'. I expect it doesn't like the ">". Quoting the SQL and bracketing the where clause portion doesn't help. Perhaps there's an escape sequence I could use? And is there point in constructing a view with columns of WKT geometry, XMIN, YMIN, XMAX and YMAX to see if the ODBC driver would take advantage of those? Thanks! Brent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR and database filter on spatial extent
More performance news... Using ogrinfo -> VRT -> ODBC -> SqlServer2000, the query using: alta83_v1 - takes less than 2 seconds Select * from alta83_v1 - a very simple select, logically equivalent to the above - takes 30 seconds (!) From the debug messages below, it seems that the spatial filter is not passed to the database server if a is included in the VRT. Not too surprising given the difficultly of parsing a SELECT statement and rebuilding it with additional clauses, but we should mention the order of magnitude performance hit in the doc. And I'm not sure why ogr seems to execute the SQL three times... Thanks! Brent = alta83_v1: - G:\GeoData\Wells>ogrinfo -spat -114 51 -113 52 -ro --debug ON all_wells.ovf wells | more OGR_ODBC: EstablishSession(DSN:"Wells_SQL", userid:"test", password:"test") ODBC: SQLConnect(Wells_SQL) OGR_ODBC: Table alta83_v1 has no identified FID column. OGR: OGROpen(ODBC:test/t...@wells_sql,alta83_v1/00C846D0) succeeded as ODBC. OGR: OGROpen(all_wells.ovf/00C83B60) succeeded as VRT. OGR: GetLayerCount() = 1 OGR_ODBC: ExecuteSQL(SELECT * FROM alta83_v1 WHERE BOT_LONG > -114 AND BOT_LONG < -113 AND BOT_LAT > 51 AND BOT_LAT < 52) OGR_ODBC: ExecuteSQL(SELECT * FROM alta83_v1 WHERE BOT_LONG > -114 AND BOT_LONG < -113 AND BOT_LAT > 51 AND BOT_LAT < 52) OGR_ODBC: ExecuteSQL(SELECT * FROM alta83_v1 WHERE BOT_LONG > -114 AND BOT_LONG < -113 AND BOT_LAT > 51 AND BOT_LAT < 52) INFO: Open of `all_wells.ovf' using driver `VRT' successful. Select * from alta83_v1 G:\GeoData\Wells>ogrinfo -spat -114 51 -113 52 -ro --debug ON all_wells.ovf wells | more OGR_ODBC: EstablishSession(DSN:"Wells_SQL", userid:"test", password:"test") ODBC: SQLConnect(Wells_SQL) OGR_ODBC: Table alta83_v1 has no identified FID column. OGR: OGROpen(ODBC:test/t...@wells_sql,alta83_v1/00C846D0) succeeded as ODBC. ODBC: ExecuteSQL(Select * from alta83_v1) called. OGR_ODBC: Table SELECT has no identified FID column. OGR: OGROpen(all_wells.ovf/00C83B60) succeeded as VRT. OGR: GetLayerCount() = 1 OGR_ODBC: Recreating statement. OGR_ODBC: Recreating statement. INFO: Open of `all_wells.ovf' using driver `VRT' successful. Brent Fraser wrote: Frank, Interesting. I'm attempting to use Mapserver -> OGR -> VRT -> ODBC -> SqlServer2000 -> a table of 250k rows of point features. Testing with ogrinfo to return all the rows takes about 25 seconds, but with a "-spat" to get about 3000 rows takes < 2 seconds (very acceptable). And I didn't need to create a view with XMIN,YMIN,XMAX,YMAX as implied on GDAL's ODBC format page; OGR applied the filter to the point coordinate columns (dunno what it would do with lines or polygons and no min/max columns though). Looking at mapserver's mapogr.cpp, it appears that OGR_L_SetSpatialFilter is being called, so I'll have to do some more tracing/debugging to find out why my mapserver performance is so bad with this layer. Thanks! Brent Fraser Frank Warmerdam wrote: Brent Fraser wrote: Hi All, In the case of accessing data in a relational database, does OGR have the ability to pass a spatial extent to the database to use as a filter on the geometry before sending the rows? Brent, Yes. The OGRLayer has a SetSpatialFilter() method for this. Some drivers evaluate the spatial filter in OGR after reading all records, but smart drivers are able to use the spatial filter for efficient querying. So, the spatially enabled databases work it into the query. The -spat switch for ogrinfo is translated into a SetSpatialFilter() call for instance. Best 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] OGR and database filter on spatial extent
Frank, Interesting. I'm attempting to use Mapserver -> OGR -> VRT -> ODBC -> SqlServer2000 -> a table of 250k rows of point features. Testing with ogrinfo to return all the rows takes about 25 seconds, but with a "-spat" to get about 3000 rows takes < 2 seconds (very acceptable). And I didn't need to create a view with XMIN,YMIN,XMAX,YMAX as implied on GDAL's ODBC format page; OGR applied the filter to the point coordinate columns (dunno what it would do with lines or polygons and no min/max columns though). Looking at mapserver's mapogr.cpp, it appears that OGR_L_SetSpatialFilter is being called, so I'll have to do some more tracing/debugging to find out why my mapserver performance is so bad with this layer. Thanks! Brent Fraser Frank Warmerdam wrote: Brent Fraser wrote: Hi All, In the case of accessing data in a relational database, does OGR have the ability to pass a spatial extent to the database to use as a filter on the geometry before sending the rows? Brent, Yes. The OGRLayer has a SetSpatialFilter() method for this. Some drivers evaluate the spatial filter in OGR after reading all records, but smart drivers are able to use the spatial filter for efficient querying. So, the spatially enabled databases work it into the query. The -spat switch for ogrinfo is translated into a SetSpatialFilter() call for instance. Best regards, ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] OGR and database filter on spatial extent
Hi All, In the case of accessing data in a relational database, does OGR have the ability to pass a spatial extent to the database to use as a filter on the geometry before sending the rows? Thanks, Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Too much precision in WKT
Even, I'll likely do the hack for now. I need to write an easy-to-use app that spits out the first and last point coordinates of the linestrings as well as the linestrings in WKT for eventual loading into a database, so I might as well hard code the precision for now. Personally, I think something like an option of: -WKT_cp eg: -WKT_cp 7 would result in 53.1467912 to specify the WKT coordinate precision would be good enough for a long term solution. Easy to implement, easy to document, and easy to check the user input. And would likely handle most of the practical precision needs. Best Regards, Brent Fraser Even Rouault wrote: Brent, no, there's currently no way to limit the precision. By looking at your example, it seems that the extra figures are significant (but perhaps not for your use case). You'd get 01 or 99 at the end of the numbers if they were not significant So there's no way OGR can guess that you don't want them. The WKT formatting for a point is done by OGRMakeWktCoordinate() in ogr/ogrutils.cpp and it always outputs 15 digits after the comma. This is the place where to hack if you want to change that default behaviour. The easiest way to add this would be adding a configuration option. The other - and cleaner - possibility is to propagate a precision parameter into the whole sequence of calls of ExportToWkt(). In any case, care should be taken that the output buffer passed by calling functions is sufficiently large. Currently all callers are supposed to pass a sufficiently large buffer, but if precision is made parametrable, this would be a bit more complicated : that wouldn't be a bad thing by the way to revisit that, as currently the code that allocates the buffer has plenty of magic - and different - constants all over and different tests to check if it must be resized... But specifying the formatting is not as trivial as one could think at first. For a floating point, you have several possibilities that are reflected by the possibilities offered by printf with floating numbers : - fixed point notiation (%f) where you can control the maximum number of figures after the comma, but not before - exponential notation (%e) where you can control the number of significant figures. - fixed point or exponential (%g), whichever is more appropriate for its magnitude I'm not sure how the user could specify the format he wants. Passing the printf formatting string is a bit dangerous, as very bad things could happen if he got wrong... I found that Frank had written a paragraph somehow a bit related to that issue, but it was more about OGC WKT Coordinate System. See "Numerical Precision in WKT" http://home.gdal.org/projects/opengis/wktproblems.html I also found the following discussion : http://lists.ingres.com/pipermail/gis-users/2009-January/000145.html Best regards, Even P.S : I had a doubt if the OGC WKT spec allowed exponential notation, but apparently yes (OGC 06-103r3 page 53-54, OpenGIS® Implementation Specification for Geographic information - Simple feature access - Part 1: Common architecturen v1.2). Le Tuesday 02 June 2009 19:50:40 Brent Fraser, vous avez écrit : I've been experimenting with v1.6.0 ogr2ogr: ogr2ogr -f csv test_dir test_in.shp -nln test_out -lco GEOMETRY=AS_WKT The precision of the coordinates in the WKT seems to be overkill, eg: "LINESTRING (-115.11433812265155 53.146791166875367,-115.12192424362472 53.147304268559473, Is there any way to limit the precision? Thanks! Brent Fraser ___ 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] Too much precision in WKT
Frank, I'm all for that approach. The option should be added to ogrinfo as well since it can dump out WKT when incanted properly (I had started to sed/grep/awk its output, then I found that ogr2ogr -f CSV would dump WKT geometry). Many Thanks, Brent Fraser Frank Warmerdam wrote: Brent Fraser wrote: I've been experimenting with v1.6.0 ogr2ogr: ogr2ogr -f csv test_dir test_in.shp -nln test_out -lco GEOMETRY=AS_WKT The precision of the coordinates in the WKT seems to be overkill, eg: "LINESTRING (-115.11433812265155 53.146791166875367,-115.12192424362472 53.147304268559473, Is there any way to limit the precision? Brent, This question has come up a few times in the context of different text oriented drivers. I'm wondering if we should add an option to ogr2ogr (and a method on OGRGeometry) to reduce coordinate precision to a specified amount. I think I prefer this approach to the alternative of adding creation options or configuration variables to each driver text driver. Best regards, ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Too much precision in WKT
Peter, I want to round the coordinates of the vertices of linestrings represented as WKT, not the projection parameters. My problem is that databases have different field types for strings depending on their expected length (MS Access TEXT vs MEMO for example, and I'm still trying to find info on DBF types). As you've explained below, instead of representing the latitude of a vertex as 53.146791166875367, it would be ok to instead use 53.1467912 and be within a centimeter of the original latitude (as long as the projection parameters remained at their full precision). Best Regards, Brent Fraser Peter J Halls wrote: Brent, why would you want to? Maybe you do not appreciate the implications of so doing? The parameters of which you complain define the ellipsoid, the shape and measurements of the Earth, to be used for that projection. Round them and you effectively change the size and shape of the Earth! Based on the Clarke Ellipsoid, 1m at the equator is approximately 0.0008833141949 of a degree; round that and you lose the positional accuracy - and 1m positional accuracy is, frankly, not great. You can divide that figure by ten if you want to work at .1m positional accuracy. However, most mapping agencies and surveyors will probably use 1cm positional accuracy, or one hundreth: 0.08833141949. There are a number of other spheroids and ellipsoids in use, as the surface of our planet is sufficiently irregular to necessitate different measurements depending upon the location of interest. The map projection equations are sensitive to this numerical precision - they have to be - so reducing the numerical precision of the projection parameters will significantly reduce the precision of the final result. Indeed, it does not take much rounding to have a dramatic result, should you seek to align the results with material projected using the proper parameters. You can test this yourself be experimenting with some spherical trigonometry - the mathematics necessary for working with position on a sphere. You can find the equations for the various projections in J P Snyder's 'Map Projections: a working manual', which is available online from USGS. Try out some of the equations by hand, with the full precision and with your proposed precision for some locations in your area of interest and I'm sure you will then appreciate the need for the precision in these parameters. Best wishes, Peter Brent Fraser wrote: I've been experimenting with v1.6.0 ogr2ogr: ogr2ogr -f csv test_dir test_in.shp -nln test_out -lco GEOMETRY=AS_WKT The precision of the coordinates in the WKT seems to be overkill, eg: "LINESTRING (-115.11433812265155 53.146791166875367,-115.12192424362472 53.147304268559473, Is there any way to limit the precision? Thanks! Brent Fraser ___ 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] Too much precision in WKT
I've been experimenting with v1.6.0 ogr2ogr: ogr2ogr -f csv test_dir test_in.shp -nln test_out -lco GEOMETRY=AS_WKT The precision of the coordinates in the WKT seems to be overkill, eg: "LINESTRING (-115.11433812265155 53.146791166875367,-115.12192424362472 53.147304268559473, Is there any way to limit the precision? Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] more commands, less utilities
How about gdaldem2image for a name? A little long, but less confusion over what it actually does... Matt Wilkie wrote: Upon reflection, I think we should congregate these utilities into a single utility -- gdaldem (I'm totally flexible on the name -- and have each be a separate operation within that rather than proliferating five more GDAL utilities. -- http://trac.osgeo.org/gdal/ticket/2640 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Splitting and projecting a large .sid image
>> 5 Use gdal_translate with the extents from step 2 to extract the output >> tile pixels from its corresponding vrt. > > H, where do you reproject the raster images? I thought you had to > use gdalwarp for that. Or is that step 6? Oops. Yep, use gdalwarp (I was looking at my DEM processing script; its geo-to-geo). Brent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Splitting and projecting a large .sid image
Steve, My general strategy when doing this (on terabytes of raster data) is: 1. Create a tileindex of input data. Since you've got multiple zones of raster data, you'll need to create one tileindex for each zone (using gdaltindex), de-project them to geographic, and concatenate them (use ogr2ogr or PostGIS). 2. Create a geographic tileindex of output files (via your own script) The size, extents, and amount of overlap are up to you. It doesn't have to be a tileindex; it could be just a list of extents for use in steps 3 and 5 below. 3. Do a spatial select for each output tile on the input tileindex to create a "mini" tileindex for each output tile. I use a hacked version of shputils to do this, but I expect PostGIS or ogr2ogr + SQL would be a better way to go. 4. Convert all the mini tileindexes to vrt files (see the new gdalbuildvrt app!) 5 Use gdal_translate with the extents from step 2 to extract the output tile pixels from its corresponding vrt. Convoluted yes, but very "scalable". And no Null pixels in the output tiles. Brent Fraser GeoAnalytic Inc. Stephen Woodbridge wrote: Hi all, This is my annual wake-up got to do something with rasters and forgot the finer details of the last time I did a problem similar to it. I have a bunch of LandSat mrsid imagery, they are in multiple UTM zone projections and I want to convert them to geographic projection, WGS84, and into tiff files. The files are too large to convert into a single tif, so I plan to quarter each file into 4 pieces. So I wrote a script the reads the files and collects info about them, then generates 4 gdalwarp commands to reproject each quarter of the file (perl included below). A typical command looks like: gdalwarp -srcnodata 0 -dstnodata 0 -s_srs +init=epsg:32642 -t_srs EPSG:4326 -te 68.676729652 24.9707453133203 72.7715338115736 27.4799559407707 -rb -wm 250 --config GDAL_ONE_BIG_READ ON -co "TILED=YES" ./GeoCover2000/n-42-25.sid ./af/n-42-25_1-0.tif The reprojection seems to have worked fine, but the bounds seem to be off because of the black (nodata) areas. You can see my poor results here: http://imaptools.com/downloads/raster.gif Is there a better way to do this? I seem to remember a tool to generate tiles. Does this need to be done in two steps? How? Did I mess up my logic below? I've been over it multiple times and do not see any errors. Any help would be appreciated. Thanks, -Steve $ gdalinfo -nomd ./GeoCover2000/n-42-25.sid Driver: MrSID/Multi-resolution Seamless Image Database (MrSID) Files: ./GeoCover2000/n-42-25.sid Size is 51078, 39090 Coordinate System is `' Origin = (136066.125,3323577.375) Pixel Size = (14.250,-14.250) Corner Coordinates: Upper Left ( 136066.125, 3323577.375) Lower Left ( 136066.125, 2766544.875) Upper Right ( 863927.625, 3323577.375) Lower Right ( 863927.625, 2766544.875) Center ( 46.875, 3045061.125) Band 1 Block=1024x128 Type=Byte, ColorInterp=Red Minimum=0.000, Maximum=226.000, Mean=109.610, StdDev=39.202 Overviews: 25539x19545, 12770x9773, 6385x4887, 3193x2444, 1597x1222, 799x611, 400x306, 200x153 Band 2 Block=1024x128 Type=Byte, ColorInterp=Green Minimum=0.000, Maximum=229.000, Mean=125.833, StdDev=31.828 Overviews: 25539x19545, 12770x9773, 6385x4887, 3193x2444, 1597x1222, 799x611, 400x306, 200x153 Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue Minimum=0.000, Maximum=227.000, Mean=107.092, StdDev=37.230 Overviews: 25539x19545, 12770x9773, 6385x4887, 3193x2444, 1597x1222, 799x611, 400x306, 200x153 $ cat ./GeoCover2000/n-42-25.prj Projection UTM Datum WGS84 Zunits NO Units METERS Zone 42 Xshift 0.00 Yshift 0.00 Parameters Perl snippet: foreach my $sid (@files) { next unless $wanted && $sid =~ /$wanted/i; if (! -r $sid) { warn "SKIPPING: $sid is not readable!\n"; next; } my $gdalinfo = `gdalinfo -nomd $sid`; # parse out the strings that we need $gdalinfo =~ /Upper Left \((\s*[^,]+),\s*([^\)]+)/s || warn "Upper Left parse error\n"; my $ul = [$1, $2]; $gdalinfo =~ /Lower Left \((\s*[^,]+),\s*([^\)]+)/s || warn "Lower Left parse error\n"; my $ll = [$1, $2]; $gdalinfo =~ /Upper Right \((\s*[^,]+),\s*([^\)]+)/s || warn "Upper Right error\n"; my $ur = [$1, $2]; $gdalinfo =~ /Lower Right \((\s*[^,]+),\s*([^\)]+)/s || warn "Lower Right error\n"; my $lr = [$1, $2]; my $outf = $sid; $outf =~ s/\.sid//; $outf =~ m/\/([^\/]+)$/; $outf = $1; my $prj = $sid; $prj =~ s/\.sid/.prj/; if (! -r $prj) { warn "SKIPPING: $prj is not readable!\n"; next; }
Re: [gdal-dev] Transparency using gdalbuildvrt
Jason, I'm doing some similar things with vrt and DEMs. I can't comment about the alpha/RGBA thing but I was able to get the no-data setting used for the source data by editing the vrt file. gdalbuildvrt creates a "SimpleSource" type of VRT: - -3.276700E+004 082\082a13_0102_deme.dem 1 : - I edited the VRT to be a "ComplexSource" type, and added a tag for each source entry: - -3.276700E+004 082\082a13_0102_deme.dem 1 -32767 : ----- Best Regards, Brent Fraser Jason Beverage wrote: Hi all, I've been playing around with gdalbuildvrt lately and it has proven to be a very useful tool:) What I'm attempting to do is take multiple RGB images and treat them as a single image using the VRT. This works great, but since the images are RGB, areas where there is no data are coming back as black. Is there a way that I can modify the VRT or the code to treat the VRT as an RGBA dataset and somehow say "If there is no data in any bands, the alpha is 0"? Thanks! Jason ___ 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] gdalbuildvrt and location, Location, LOCATION
Yep, my fault. gdaltindex and gdalbuiltvrt are just fine the way they are... Sorry for the noise, Brent Brent Fraser wrote: Hmm, my gdaltindex creates a DBF with "LOCATION", even though I see that gdaltindex.c sets the default to "location" on line 71: const char *tile_index = "location"; I'll have to step through the code with the debugger... Thanks! Brent Even Rouault wrote: By default, gdaltindex produces a DBF with 'location' (lowercase) and ogrtindex 'LOCATION' (uppercase). So I think I've added this check so that people don't try to use ogrtindex output. The '-tileindex' option is indeed done for those specific cases where non-standard GDAL tileindex must be used. Le Thursday 12 March 2009 19:58:14 Brent Fraser, vous avez écrit : So either the message should be changed (since the index WAS built using GDAL products), or the default should be changed to LOCATION (my preference). I'd be happy to file a bug report if someone has a strong feeling either way.. ___ 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] gdalbuildvrt and location, Location, LOCATION
Hmm, my gdaltindex creates a DBF with "LOCATION", even though I see that gdaltindex.c sets the default to "location" on line 71: const char *tile_index = "location"; I'll have to step through the code with the debugger... Thanks! Brent Even Rouault wrote: By default, gdaltindex produces a DBF with 'location' (lowercase) and ogrtindex 'LOCATION' (uppercase). So I think I've added this check so that people don't try to use ogrtindex output. The '-tileindex' option is indeed done for those specific cases where non-standard GDAL tileindex must be used. Le Thursday 12 March 2009 19:58:14 Brent Fraser, vous avez écrit : So either the message should be changed (since the index WAS built using GDAL products), or the default should be changed to LOCATION (my preference). I'd be happy to file a bug report if someone has a strong feeling either way.. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] gdalbuildvrt and location, Location, LOCATION
I'm trying to use gdalbuildvrt (v1.6.0) to create a vrt file from a tileindex shapefile. The tileindex shapefile was created using gdaltindex, and dbfinfo reports: 1 Columns, 16 Records in file LOCATION string (254,0) Doing: gdalbuildvrt 082h04w2.vrt 082h04w.shp reports: This shapefile seems to be a tile index of OGR features and not GDAL products. Unable to find field `location' in DBF file `082h04w.shp'. But when doing: gdalbuildvrt -tileindex LOCATION 082h04w2.vrt 082h04w.shp I successfully get my vrt file. So either the message should be changed (since the index WAS built using GDAL products), or the default should be changed to LOCATION (my preference). I'd be happy to file a bug report if someone has a strong feeling either way.. Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Jan, I stepped through your procedure with a Canadian topo map to slightly adjust its positioning, using GDAL 1.6.0: 1. Warp raster image of topo map to a projected coordinate system No Need. Canadian topos already are georeferenced so I simply downloaded http://ftp2.cits.rncan.gc.ca/pub/canmatrix/50k_300dpi/082/h/canmatrix_082h04_tif.zip 2. Attach my new GCPs to the original topo map: gdal_translate --optfile gcps.opt -a_srs "EPSG:26912" 082h04_1_1.tif 082h04_gcps.tif (the GCPs are in UTM zone 12, as indicated by the -a_srs parameter) 3. Warp/resample image using GCPs gdalwarp -tps 082h04_gcps.tif 082h04_warped.tif 4. Examine SRS of output image: gdalinfo 082h04_warped.tif Driver: GTiff/GeoTIFF Files: 082h04_warped.tif Size is 11786, 8881 Coordinate System is: PROJCS["NAD83 / UTM zone 12N", GEOGCS["NAD83", DATUM["North_American_Datum_1983", SPHEROID["GRS 1980",6378137,298.2572221010002, AUTHORITY["EPSG","7019"]], AUTHORITY["EPSG","6269"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4269"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-111], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",50], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","26912"]] Origin = (275804.27707260562,5462603.41083706920) Pixel Size = (4.233604626074191,-4.233604626074191) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 275804.277, 5462603.411) (114d 4'56.70"W, 49d16'30.13"N) Lower Left ( 275804.277, 5425004.768) (114d 3'41.57"W, 48d56'14.30"N) Upper Right ( 325701.541, 5462603.411) (113d23'49.60"W, 49d17'28.67"N) Lower Right ( 325701.541, 5425004.768) (113d22'51.13"W, 48d57'12.15"N) Center ( 300752.909, 5443804.089) (113d43'49.87"W, 49d 6'53.14"N) Band 1 Block=11786x1 Type=Byte, ColorInterp=Palette Color Table (RGB with 256 entries) The result of step 2 (gdal_translate) DOES remove the image's original georeferencing. Actually, it replaces it with a set of GCPs (and their SRS). To leave the original would resulting in two sets of conflicting georeferencing information. And to me it doesn't matter anyway, as that "GCP-ed" image is just temporary. After the gdalwarp step I get an image with more traditional/well-supported georeferencing, I can delete it. Should gdalinfo (and possibly other programs) use the GCPs to calculate the image extents? Maybe. I'd find it useful as a quality-control check to make sure I didn't mis-code a GCP resulting in a huge raster after warping. Do you have another reason? Brent Jan Hartmann wrote: Yes exactly! Brent highlights a point I never mentioned, and that may be the source of the confusion: when I add control points to an already georeferenced image, I gdalwarp it with the -tps option: a thin spline transformer, that transforms the control points to the exact georeferenced position, but stretches the areas in between. AFAIK, this is done by triangulation. This is different from a traditional warp, in which the eventual position of the control points on the georeferenced map has errors, although minimized ones. So in matching roads on two adjacent images, the road edges on both images are assigned the same georeferenced coordinates and then, after gdalwarp -tps, they exactly align. So my workflow is: gdalwarp raw image --> georeferenced image with correct projection, based on the corner points of the original map gdal_translate : add control points to projected image, essentially small corrections of road edges and triangulation points for which we know the exact coordinates from later sources gdalwarp -tps georeferenced image --> georeferenced image with small corrections, where roads, churcht towers are on their exact locations The whole procedure is something I do regularly with ArcGIS, not only with historical maps but also with things like aerial photographs, just like Brent describes. It can be done with gdal too, as described above, if only the second step would work. I still don't see why adding control points to an image that already has a projection defined, should destroy that projection information, including its georeferenced boundaries and dimensions. Even when this is not the case: when a -a_srs flag is allowed for gdal_translate, it should be possible to add the complete geolocation informa
Re: [gdal-dev] Re: [Fwd: Re: Warping an already georeferenced image with control points
Jan, I think what you want gdalwarp to do is [Delauney] image triangulation. Not an uncommon (or unreasonable) request. I occasionally have a need for triangulation when mosaicking several satellite images. Its visually pleasing to have the roads match up precisely where the images overlap... Brent Fraser Even Rouault wrote: Jan, (I'm CC'ing the list) I'm not sure what you mean with adapting the pixelsize, now that the output has only GCPs and no more geotransform matrix. As far as including this option in baseline gdal_translate.cpp, I'm currently not really enthousiastic, as there are lots of way to get the result achieved (what is done with that patch could be done similarly with GDAL python API for instance) and the choice of corners is something rather arbitrary. Why should you trust the georeferenced coordinates of the corner points as reliable GCPs ? Why not some regularly discretized points along the edges ? Or a regular grid over the whole raster ? etc. etc. So, I'd prefer hearing from other GDAL developers or users that this patch is a sensical addition before commiting it. Even Le Thursday 12 February 2009 11:47:00, vous avez écrit : Hi Even, Thanks a lo. Just one small addition: the pixelsize has also to be adapted, else you get an image with false dimension. If it's not too much to add that, I'm going to try it out and will let you know the results. If it works, I have indeed a practical solution for my problem. Even so, I am not convinced that this is an uncommon problem. Consider a very basic GIS functionality: edge matching. You have digitized a map from many sheets, and at the end the georeferenced images don't match precisely at the edges. Edge-matching functions (ArcGIS has lots of them) ensure that the edge of a feature on one map perfectly matches the feature edge on the adjacent map. This is exactly the kind of thing you can do with gdalwarp using control points on a georeferenced image. It is also more or less the thing I want to do with my historical maps. Would this be an argument to preserve the boundaries of a georeferenced image, when adding control points? Jan Even Rouault wrote: Jan, I concur with Jukka's analysis. Your need is rather uncommon, and I also thinks it shouldn't be a default behaviour. However, adding that capability to gdal_translate is quite easy. I've attached a patch for gdal_translate.cpp that adds the capability of adding the 4 corners as 4 additional GCPs with the "-add_corners_as_gcps" option. You might give it at try (it applies cleanly on GDAL 1.6.0 and trunk). I'm not sure yet if it is valuable enough to include it in baseline gdal_translate.cpp. Best regards Le Wednesday 11 February 2009 15:59:48 Jukka Rahkonen, vous avez écrit : Jan Hartmann uva.nl> writes: Sorry to keep moaning about this, but I need an indication what's going on here. Mind, I don't need an immediate solution: for the time being I have a workaround. Just an idea whether this a real problem, a dumb question, something that can be handled in the foreseeable future (or not), perhaps with adequate funding. Everything is better than talking to a blind wall. Sorry again Frank, Jan Hi, I think that your need to keep the original extents but add more ground control points inside the image frame is rather uncommon. Doesn't it mean that you trust that the image corners are correctly warped to a new projection, but there are local distortions in the middle of image which should be corrected with a few extra ground control points? For my mind it shoudn't be the default behaviour of gdal but it might be usable as an user selectable option sometimes. I know that missing gcp's at the image corners often leads to very odd result with polynomial warping because the formula shoots over. Even unaccurately placed gcp's could help a lot in preserving the original shape of the map. I guess that you are perhaps playing with scanned historical maps? I have a few old scanned parcel maps which are covering just the area of the farm, and two of the map corners are just white background. It is impossible to measure any real ground control points from the corners because there is nothing on the map to compare with, and warping with all gcp's on the middle area makes really funny looking curves into the hand drawn rectangle framing the original map. -Jukka- ___ 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] Zonal statistics with GDAL/StarSpan
Is the FID > supposed to be in the .dbf file that accompanies the .shp file? > I have examined this file, and I don't see an FID column (column > names are in the first row). In particular, since starspan2 > reports an FID of 0 for Zemmour, Mauritania, I would expect to > see a column with the value 0 in the Zemmour row of the .dbf > file. But, there isn't one. I was assuming that FIDs were in > one of the other files in the set (.dbf, .fix, .qix, .shp, .shx). I think StarSpan will "assign" a Feature ID (FID), likely sequentially with the first one being 0, to the records as it reads them. I'm not sure if it would use an FID field in the DBF file (where are the attributes are) if it found one. Brent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Zonal statistics with GDAL/StarSpan
> Hi, > > I performed the OpenJump procedure. But, I still get: > > Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none > --vector afadmn2n.shp --fid 3 --raster rfe_2006_04_pct.tif --stats > results.txt avg --out-prefix bar --out-type table > Number of features: 550 > starspan_csv: 1: Extracting from rfe_2006_04_pct.tif > 0% terminate called after throwing an instance of > 'geos::util::TopologyException' > what(): TopologyException: side location conflict -8.25 23.1536 > Abort trap > So: > > a) Is there some way to get a dump of FIDs with their corresponding > country/province names? You could use dbfdump. I think it is delivered with FWTools. > b) What other steps might I take to detect/fix any side location > conflict(s) I've attached the source to a GDAL application I wrote (it simply accesses the GDAL/Geos buffer function) to repair geometry of polygons. Use the -r option (the -c option doesn't work). You'll need to link with a version of GDAL that uses GEOS. > c) The docs on the starspan wiki are out of sync with the current > version (1.2.03). Are there updated docs for starspan? Should I fall > back to a previous version? > I used version 1.2.01, but I don't recall why... Brent/* * SOURCE HEADER: ** Name: Purpose: Invocation: Arguments: Data Files Accessed: Program History: v1.019 Jun 2007BWF initial coding Note: Algorithm Description: Future Enhancements: *** End of SOURCE HEADER ** */ #include #include "GDAL.h" #include "cpl_string.h" #include "ogrsf_frmts.h" // // // static void Usage() { printf( "ogrClean v1.0\n" "Create clean planar polygon layer from overlapping polygons.\n" "Usage: ogrClean -a Assigned unique integer to FID attribute\n" " ogrClean -r Repair topology\n" " ogrClean -c Create non-overlapping topology\n" " ogrClean -d Delete records with NULL geometry.\n" ); exit( 1 ); } // // void notice(const char *fmt, ...) { va_list ap; fprintf( stdout, "NOTICE: "); va_start (ap, fmt); vfprintf( stdout, fmt, ap); va_end(ap); fprintf( stdout, "\n" ); } void log_and_exit(const char *fmt, ...) { va_list ap; fprintf( stdout, "ERROR: "); va_start (ap, fmt); vfprintf( stdout, fmt, ap); va_end(ap); fprintf( stdout, "\n" ); exit(1); } // // int repairPolygons(const char *pszSrcFilename, bool bOpAssignFID, bool bOpRepairTopo ) { OGRDataSource *poDS; OGRLayer*poLayer; OGRFeature *poFeature; OGRGeometry *poGeometry; OGRwkbGeometryType poGeomType; char *pszGeom=NULL; OGRErr oErr; long nLayers=0; long nFeatures=0; bool bFoundFID; int iField; long i=0, id=0; long iCountFixed=0; bool bGeomUpdated = false; /*---*/ /* Open input: */ /*---*/ OGRRegisterAll(); poDS = OGRSFDriverRegistrar::Open( pszSrcFilename, TRUE ); // we will be updating the dataset if( poDS == NULL ) { return(1); //can't open dataset } nLayers = poDS->GetLayerCount(); if (nLayers == 0) { return(2); //no layers in dataset } // TBD: allow mult layers: poLayer = poDS->GetLayer( 0 ); poLayer->ResetReading(); nFeatures = poLayer->GetFeatureCount(); if ( (poFeature = poLayer->GetNextFeature()) != NULL ) { //--- Check Geometry type: ---// poGeometry = poFeature->GetGeometryRef(); if( poGeometry != NULL ) { poGeomType = poGeometry->getGeometryType(); // TBD: check if polgon if ( poGeomType != wkbMultiPolygon && poGeomType != wkbPolygon ) { printf("Data must be polygons.\n"); } }else { return(4); //no geometry
Re: [gdal-dev] Zonal statistics with GDAL/StarSpan
Greg, OpenJump is a good topology/geometry fixing tool. Load your shapefile and do Tools->Generate->Buffer, set the Buffer distance to 0 to clean up the polygon vertices. I seem to recall it may have a problem preserving attributes when the results are saved as a shapefile, but you may be able to just use the new .shp portion with the old .dbf and .shx if the resulting geometry records are not re-ordered. Brent Greg Ederer wrote: Brent, My Google search didn't turn up much of use (most of the results were PostGIS-specific). Do you know of a tool that would tell me which feature(s) have problems? Do you know of a tool that might be able to fix these polygons? Thanks! Greg Brent Fraser wrote: Greg, It means the GEOS topology library is having problems with geometry of your input shapefile. Most likely related to duplicate vertices (or near-duplicate) in a polygon. Do a Google for "side location conflict". It would be nice if it spat out the FID of the offending polygon... Brent Greg Ederer wrote: Hi Brent, Thanks for the example. Here's what I got: Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none --vector /Users/gregederer/servers/geoserver_data/data/shapefiles/Admin/afadmn2n.shp --raster rfe_2006_04_pct.tif --stats results.txt avg starspan: --out-prefix: ? Type `starspan --help' for help I then added a --out-prefix, like so: Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none --vector /Users/gregederer/servers/geoserver_data/data/shapefiles/Admin/afadmn2n.shp --raster rfe_2006_04_pct.tif --stats results.txt avg --out-prefix bar starspan: --out-type: ? Type `starspan --help' for help I then added '--out-type table' (I found in the source that valid --out-types are: table, mini_raster_strip, mini_rasters, rasterization): Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none --vector /Users/gregederer/servers/geoserver_data/data/shapefiles/Admin/afadmn2n.shp --raster rfe_2006_04_pct.tif --stats results.txt avg --out-prefix bar --out-type table Number of features: 550 starspan_csv: 1: Extracting from rfe_2006_04_pct.tif 0% terminate called after throwing an instance of 'geos::util::TopologyException' what(): TopologyException: side location conflict -8.25 23.1536 Abort trap Any idea what the TopologyException indicates? Thanks! Greg Brent Fraser wrote: Greg, I've used StarSpan to do similar things; in my case I wanted "mode" not "avg". Using v1.2.01 (I don't think I made major mods to the source): starspan2.exe --progress --RID none --vector boundaries.shp --raster big.tif --stats results.txt mode I edited the resulting results.txt text (csv) file and used GDAL's ogr2ogr to do a join on the shapefile's DBF file based on the FID. Brent Fraser Greg Ederer wrote: Hello all, I'm a GDAL newbie. I have a shape file containing subnational boundaries for Africa. I have a raster containing rainfall for Africa. I need to pull out the spatial average of rainfall for each of the subnational units. I think that I might be able to pull this off using StarSpan. But, the docs on the StarSpan wiki do not line up with the current version of the program (the changelog mentions that certain command line switches have been eliminated, and others added; but it doesn't say what the new switches do, or how to use them). I am digging through the StarSpan sources to try to figure out how to use it. Does anyone have any suggestions as to how I might get the needed data using GDAL (with, or without StarSpan)? Thanks! Greg ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Zonal statistics with GDAL/StarSpan
Greg, It means the GEOS topology library is having problems with geometry of your input shapefile. Most likely related to duplicate vertices (or near-duplicate) in a polygon. Do a Google for "side location conflict". It would be nice if it spat out the FID of the offending polygon... Brent Greg Ederer wrote: Hi Brent, Thanks for the example. Here's what I got: Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none --vector /Users/gregederer/servers/geoserver_data/data/shapefiles/Admin/afadmn2n.shp --raster rfe_2006_04_pct.tif --stats results.txt avg starspan: --out-prefix: ? Type `starspan --help' for help I then added a --out-prefix, like so: Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none --vector /Users/gregederer/servers/geoserver_data/data/shapefiles/Admin/afadmn2n.shp --raster rfe_2006_04_pct.tif --stats results.txt avg --out-prefix bar starspan: --out-type: ? Type `starspan --help' for help I then added '--out-type table' (I found in the source that valid --out-types are: table, mini_raster_strip, mini_rasters, rasterization): Macintosh-4:starspan_tests gregederer$ starspan2 --progress --RID none --vector /Users/gregederer/servers/geoserver_data/data/shapefiles/Admin/afadmn2n.shp --raster rfe_2006_04_pct.tif --stats results.txt avg --out-prefix bar --out-type table Number of features: 550 starspan_csv: 1: Extracting from rfe_2006_04_pct.tif 0% terminate called after throwing an instance of 'geos::util::TopologyException' what(): TopologyException: side location conflict -8.25 23.1536 Abort trap Any idea what the TopologyException indicates? Thanks! Greg Brent Fraser wrote: Greg, I've used StarSpan to do similar things; in my case I wanted "mode" not "avg". Using v1.2.01 (I don't think I made major mods to the source): starspan2.exe --progress --RID none --vector boundaries.shp --raster big.tif --stats results.txt mode I edited the resulting results.txt text (csv) file and used GDAL's ogr2ogr to do a join on the shapefile's DBF file based on the FID. Brent Fraser Greg Ederer wrote: Hello all, I'm a GDAL newbie. I have a shape file containing subnational boundaries for Africa. I have a raster containing rainfall for Africa. I need to pull out the spatial average of rainfall for each of the subnational units. I think that I might be able to pull this off using StarSpan. But, the docs on the StarSpan wiki do not line up with the current version of the program (the changelog mentions that certain command line switches have been eliminated, and others added; but it doesn't say what the new switches do, or how to use them). I am digging through the StarSpan sources to try to figure out how to use it. Does anyone have any suggestions as to how I might get the needed data using GDAL (with, or without StarSpan)? Thanks! Greg ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Zonal statistics with GDAL/StarSpan
Greg, I've used StarSpan to do similar things; in my case I wanted "mode" not "avg". Using v1.2.01 (I don't think I made major mods to the source): starspan2.exe --progress --RID none --vector boundaries.shp --raster big.tif --stats results.txt mode I edited the resulting results.txt text (csv) file and used GDAL's ogr2ogr to do a join on the shapefile's DBF file based on the FID. Brent Fraser Greg Ederer wrote: Hello all, I'm a GDAL newbie. I have a shape file containing subnational boundaries for Africa. I have a raster containing rainfall for Africa. I need to pull out the spatial average of rainfall for each of the subnational units. I think that I might be able to pull this off using StarSpan. But, the docs on the StarSpan wiki do not line up with the current version of the program (the changelog mentions that certain command line switches have been eliminated, and others added; but it doesn't say what the new switches do, or how to use them). I am digging through the StarSpan sources to try to figure out how to use it. Does anyone have any suggestions as to how I might get the needed data using GDAL (with, or without StarSpan)? Thanks! Greg ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] resampling artifacts
Even, It is fixed in v1.6 Beta. Thanks! Brent Even Rouault wrote: Brent, Have you tried with GDAL 1.6.0 beta ? Much work has been done in that area and those artifacts should have disappeared and people. See http://trac.osgeo.org/gdal/ticket/2627 Le Friday 07 November 2008 00:08:51 Brent Fraser, vous avez écrit : Has anyone experienced resampling artifacts using gdalwarp? I'm using v1.5.2, and with lanczos or cubicspline resampling, I get a vertical line and a horizontal line of seemingly random pixels about 6 pixels wide running through the center of the output image. I'm resampling a 16-bit tif image from 12.5 meters to 200 meters. The other resampling methods result in an ok output image. I've attached a png to illustrate the effect. I thought I'd ask before filing a bug report... Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] resampling artifacts
Thanks Even, I'll try the 1.6.0 beta. Even Rouault wrote: Brent, Have you tried with GDAL 1.6.0 beta ? Much work has been done in that area and those artifacts should have disappeared and people. See http://trac.osgeo.org/gdal/ticket/2627 Le Friday 07 November 2008 00:08:51 Brent Fraser, vous avez écrit : Has anyone experienced resampling artifacts using gdalwarp? I'm using v1.5.2, and with lanczos or cubicspline resampling, I get a vertical line and a horizontal line of seemingly random pixels about 6 pixels wide running through the center of the output image. I'm resampling a 16-bit tif image from 12.5 meters to 200 meters. The other resampling methods result in an ok output image. I've attached a png to illustrate the effect. I thought I'd ask before filing a bug report... Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] resampling artifacts
Has anyone experienced resampling artifacts using gdalwarp? I'm using v1.5.2, and with lanczos or cubicspline resampling, I get a vertical line and a horizontal line of seemingly random pixels about 6 pixels wide running through the center of the output image. I'm resampling a 16-bit tif image from 12.5 meters to 200 meters. The other resampling methods result in an ok output image. I've attached a png to illustrate the effect. I thought I'd ask before filing a bug report... Thanks! Brent Fraser <>___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [Gdal-dev] OGR: DXF driver
Frank Warmerdam wrote: I would note that the activity is not terribly complicated, and I would welcome others taking a crack at it. I am not particularly trying to drum up business for myself here. I should also note that there is a contributed "write only" dxf driver sitting in Trac currently waiting for integration. It is using dxflib (which I'm not so keen on), and doesn't do the read side that is of greatest interest. But it might be a base on which someone could build the read side of things. http://trac.osgeo.org/gdal/ticket/2555 Best regards, Frank, FYI, there's a minimalist DXF reader in VTP (VTP\TerrainSDK\vtdata\DxfParser.cpp). Likely a little too c++ for GDAL/OGR to be used as-is, but might be a useful supplement to the doc you've gathered (and the licensing is as libre as it gets). Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] copy GCPs from .xml to GTiff
Shawn, I had a similar need some time ago. I hacked gdalinfo to output only the GCPs in a GDAL-friendly format (one "-GCP Pixel Line X Y" per line) and redirected the results to a text file of command-line options: my_gdalinfo file1.tif > file1_gcps.opt Then attached them to my other tiff using stock gdal_translate: gdal_translate --optfile file1_gcps.opt file2.tif file2b.tif Brent Gong, Shawn (Contractor) wrote: Hi list, I have two Radarsat-2 files File A: product.xml File B: single channel GTiff I am looking for a simple way to copy file A’s GCPs and GCP Projection to file B, delete 4 existing GCPs in file B, and resave (as GTiff). Basically replacing the below Blue text with Red text. It will be even better if Pink text can also be copied over. The closest code example I found is to create a VRT wrapper around the dataset: self.src_ds = gdal.Open(‘product.xml’) opts = vrtutils.VRTCreationOptions(self.src_ds.RasterCount) lines = gdal.SerializeXMLTree(vrtutils.serializeDataset(self.src_ds, opts)) vrtstr = string.join(lines,'') vrtds = gdal.Open(vrtstr) self.new_ds = gdal.GetDriverByName("GTiff").CreateCopy(tmp_filename, vrtds, 0) read in array... gdalnumeric.BandWriteArray(self.new_ds.GetRasterBand(1), datablock, x_size, y_size) However the new file has the same number of bands as product.xml (more than one band). I also followed OpenEV's compose.py and got all the GCPs from file A into a list. But don't know how to save them to file B. Your help is appreciated. Shawn *Gdalinfo file A:* Driver: RS2/RadarSat 2 XML Product Files: product.xml Size is 12132, 3324 Coordinate System is `' GCP Projection = GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]] GCP[ 0]: Id=1, Info= (0,0) -> (-6.44497337453028,36.2346972769233,41.6078605651855) GCP[ 1]: Id=2, Info= (1213.09998,0) -> (-6.27897074066121,36.2105053710372,41.6078605651855) GCP[ 2]: Id=3, Info= (2426.19995,0) -> (-6.11306984793545,36.1860836654035,41.6078605651855) GCP[ 3]: Id=4, Info= (3639.30005,0) -> (-5.9472716724433,36.1614325255025,41.6078605651855) ... ... GCP[118]: Id=119, Info= (9704.7998,3323) -> (-5.20826990376107,35.666853278978,41.6078605651855) GCP[119]: Id=120, Info= (10917.9004,3323) -> (-5.04390384909976,35.6408291207172,41.6078605651855) GCP[120]: Id=121, Info= (12131,3323) -> (-4.87964636502838,35.614581497011,41.6078605651855) Metadata: PRODUCT_TYPE=SGF PIXEL_SPACING=1.2500e+01 LINE_SPACING=1.2500e+01 BETA_NOUGHT_LUT=lutBeta.xml SIGMA_NOUGHT_LUT=lutSigma.xml GAMMA_LUT=lutGamma.xml SATELLITE_IDENTIFIER=RADARSAT-2 SENSOR_IDENTIFIER=SAR BEAM_MODE=W2 ACQUISITION_START_TIME=2008-02-13T06:26:25.356741Z NEAR_RANGE_INCIDENCE_ANGLE=3.06394997e+01 FAR_RANGE_INCIDENCE_ANGLE=3.94798698e+01 SLANT_RANGE_NEAR_EDGE=9.078054910011414e+05 Subdatasets: SUBDATASET_3_NAME=RADARSAT_2_CALIB:BETA0:product.xml SUBDATASET_3_DESC=Beta Nought calibrated SUBDATASET_2_NAME=RADARSAT_2_CALIB:SIGMA0:product.xml SUBDATASET_2_DESC=Sigma Nought calibrated SUBDATASET_4_NAME=RADARSAT_2_CALIB:GAMMA:product.xml SUBDATASET_4_DESC=Gamma calibrated SUBDATASET_1_NAME=RADARSAT_2_CALIB:UNCALIB:product.xml SUBDATASET_1_DESC=Uncalibrated digital numbers Corner Coordinates: Upper Left (0.0,0.0) Lower Left (0.0, 3324.0) Upper Right (12132.0,0.0) Lower Right (12132.0, 3324.0) Center ( 6066.0, 1662.0) Band 1 Block=12132x43 Type=UInt16, ColorInterp=Undefined Metadata: POLARIMETRIC_INTERP=VV Band 2 Block=12132x43 Type=UInt16, ColorInterp=Undefined Metadata: POLARIMETRIC_INTERP=VH *Gdalinfo file B:* Driver: GTiff/GeoTIFF Files: imagery_VH.tif Size is 12132, 3324 Coordinate System is `' GCP Projection = GEOGCS["User-Defined",DATUM["unknown",SPHEROID["unnamed",6378137,298.257222932867]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]] GCP[ 0]: Id=1, Info= (0.5,0.5) -> (-6.44505378645208,36.234651733764,41.6078221084219) GCP[ 1]: Id=2, Info= (12131.5,0.5) -> (-4.78969807179151,35.9824333182103,41.6079517686606) GCP[ 2]: Id=3, Info= (12131.5,3323.5) -> (-4.8797286318157,35.6145371536105,41.607949202978) GCP[ 3]: Id=4, Info= (0.5,3323.5) -> (-6.52707534347629,35.8669100589582,41.6078219661592) Metadata: AREA_OR_POINT=Area TIFFTAG_IMAGEDESCRIPTION={ bandList = [ 4; ] } TIFFTAG_DATETIME=2008:02:13 22:02:16 TIFFTAG_MINSAMPLEVALUE=24 TIFFTAG_MAXSAMPLEVALUE=18827 Image S
[gdal-dev] Status of RFC 4: Geolocation Arrays especially Envisat wrt GCPs
Frank, The short story: I've got some Single Look Complex EnviSat data I need to process. Will gdalwarp do a piecewise (triangulated) warp using GCPs or only(!) use them to calculate an image-wide polynomial? The longer story: ESA's Best (http://earth.esa.int/services/best/) does an OK job extracting the amplitude from the data, but because it uses tiff files for intermediate storage, and their implementation is limited to 2gb tiffs, it's necessary to tile an input image. Best assigns a different affine transformation to each output GeoTiff tile, which causes a mis-alignment of the tiles when re-assembling them (using gdalwarp etc). Best does output an extensive metadata text file with each tile including a geolocation "array" of latitude and longitude. Unfortunately the "array" is not regular (eg, my test tile has values for lines 1, 811, 812, 1622, 1623, 2240; every 165 pixels except the last interval which is 195 pixels). If gdalwarp would use this semi-regular array to triangulate the pixel locations, that would be great, otherwise I'll have to use something else... Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] inflateCopy required?
Frank Warmerdam wrote: Brent Fraser wrote: I just did a checkout and build of GDAL, and got the following: Creating library gdal_i.lib and object gdal_i.exp cpl_vsil_gzip.obj : error LNK2019: unresolved external symbol _inflateCopy referenced in function "public: class VSIGZipHandle * __thiscall VSIGZipHandle::Duplicate(void)" ([EMAIL PROTECTED]@@[EMAIL PROTECTED]) gdal16dev.dll : fatal error LNK1120: 1 unresolved externals Where can I find the "inflateCopy " function? Zlib? Should nmake.opt be updated to include a definition to Zlib? Brent, Normally I think the internal copy of zlib (gdal/frmts/zlib) should provide that function. Is it possible you are building in an unconventional fashion or are using an external zlib of a different version? Best regards, All I did was edit nmake.opt for locations of ECW, OGDI, MrSID, HDF4, and PG libs. Hmm, there are no .obj files in the frmts\zlib dir... And the frmts\makefile.vc has: !IFDEF OGDIDIR EXTRAFLAGS = $(EXTRAFLAGS) -DFRMT_ogdi !ELSE EXTRAFLAGS = $(EXTRAFLAGS) -DFRMT_zlib !ENDIF I guess the zlib functions in my ogdi build are too old. I'll look into that... Thanks! Brent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] inflateCopy required?
I just did a checkout and build of GDAL, and got the following: Creating library gdal_i.lib and object gdal_i.exp cpl_vsil_gzip.obj : error LNK2019: unresolved external symbol _inflateCopy referenced in function "public: class VSIGZipHandle * __thiscall VSIGZipHandle::Duplicate(void)" ([EMAIL PROTECTED]@@[EMAIL PROTECTED]) gdal16dev.dll : fatal error LNK1120: 1 unresolved externals Where can I find the "inflateCopy " function? Zlib? Should nmake.opt be updated to include a definition to Zlib? Thanks! Brent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Topological Union?
Craig, I've had limited success with the Union function in Geos (via GDAL). It's slow and fails on large polygon sets (some of my files have 10,000 polygons with some convoluted geometry), but if your sets are smaller, it may work. FYI, Geos is derived (ported?) from JTS; Martin Davis has a good explanation of union strategies at http://lin-ear-th-inking.blogspot.com/2007/11/fast-polygon-merging-in-jts-using.html. Currently, I'm considering using Geotools (http://geotools.codehaus.org/) as a framework to get at the JTS functionality. And there are other polygon processing libraries out there. There's a nice summary at http://www.complex-a5.ru/polyboolean/comp.html, but I haven't tried any of those yet. Good Luck! Brent Fraser Mateusz Loskot wrote: Craig Miller wrote: Thanks Mateusz. That's exactly what I was thinking of doing if someone hasn't put it together already. Craig, BTW, you don't have to convert OGRGeometry (you play with it when reading Shapefile using OGR) to geos::geom::Geometry. You can union OGRGeometry using OGRGeometry::Union() or its C API equivalent OGR_G_Union These operations are based on GEOSUnion function. Best regards, ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Disabiling Geos messages via OGR
I'm using OGR's Geos functions (e.g. OGRGeometry::IsValid) to process some polygons. Geos spits out messages to console (stderr?) like "Warning 1: Self-intersection at or near point ...". Is there a way to disable these messages via OGR? Thanks! Brent Fraser ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev