# Tag Info

7

It is easier with Fiona, more "Pythonic", and list slicing: import fiona with fiona.drivers(): for line in fiona.open("some_shapefile.shp"): # print first and last point of every line print line['geometry']['coordinates'][0], line['geometry']['coordinates'][-1] And with shapely: from shapely.geometry import Point for line in ...

4

ogrinfo can shorten the output considerably using the -so flag. -so: Summary Only: supress listing of features, show only the summary information like projection, schema, feature count and extents. So ogrinfo -ro -so file.shp should give a summary of the metadata. And -al: List all features of all layers (used instead of having to give layer ...

3

You can do this with the GDAL/OGR python bindings. Here's a link to the OGR API tutorial. A worked example: from osgeo import ogr ds=ogr.Open(somepolylines) lyr=ds.GetLayer() for i in range(lyr.GetFeatureCount()): feat=lyr.GetFeature(i) geom=feat.GetGeometryRef() firstpoint=geom.GetPoint(0) lastpoint=geom.GetPoint(geom.GetPointCount()-1) ...

2

Based on huckfinn answers, few other comments and together with my findings: Winning format is JPEG2000 (why and which version is mentioned below Why not others) Why not others: JPEG Size limitation both data size and dimensions (4GB and 65500x65500) no (internal) pyramids possibility = bigger the image the longer it takes to display it when pan/zoom ...

2

OSGEO4W and all standalone QGIS installers come with a OSGEO4W Shell. Start that, and type gdalinfo --version and read the result. You may have different versions on the disk: Standalone, OSGEO4W and also from gisinternals if you want the latest GDAL build, but every package sets its environment so that it is using the version it was delivered with. ...

2

Try this out, if the exit status is 0 then the command has run as expected. program = "C:\Path\to\GDAL" command = 'gdal_rasterize -a ICE_TYPE -where \"ICE_TYPE=\'Open Water\'\" -b 1 -burn 5 -l ' + shapefileshortname +' ' + shapefile + ' ' + outraster input_command = [command] command_run = subprocess.call([program, input_command]) if command_run == 0: ...

2

What I would recommend, although you could also file through this with some sort of batch loop: dir /b *.tif > list_of_tiff_files.txt gdalbuildvrt list_of_tiff_files.txt vrt_of_files.vrt gdal_merge vrt_of_files.vrt merged.tif So gdalbuildvrt works with a list of files, and gdal_merge should work with a virtual raster file (vrt).

2

I recommend the use of a vrt file, a virtual raster file that would point to your data with the bounding box that you want. you can call it from python using import subprocess subprocess.call(["gdalbuildvrt", "-te", "xmin", "ymin", "xmax", "ymax", "output.vrt", "input.tif"]) then you use your vrt like an image. If you want you can create the vrt with ...

2

The simple way to do this in QGIS is to use the Raster Calculator (Raster->Raster Calculator). You have a couple of options. The easiest to explain/understand is to make a unitary raster from your mask (all data set to either 1 or NoData) and then multiply your clip layer by the unitary mask layer. To ensure the extents match the mask layer, in the ...

1

This is a one-band raster (not RGB) having pixels with values of 0 and 255, yes? Perhaps the raster would benefit from a little pre-processing, you might like to try Raster | Analysis | Sieve before using the raster polygoniser. You might also have to experiment with the sieve settings. The polygoniser doesn't know anything about colours, only about pixel ...

1

I never used QGIS for that, so I can only offer some generic advice about things to check: try starting QGIS from the command line. It may be that it reports some error or exception on the console If I am not mistaken, QGIS uses GDAL, and therefore the GDAL Georaster driver. Try setting the CPL_DEBUG environment variable to ON before running QGIS from the ...

1

Your command line will try to merge mosaic.tif (as part of *.tif) with itself. Creating the mosaic in another folder would prevent that. The correct syntax for batch processing is: for %%N in (D:\Karten\gdal\gdal2tiles\NL50\*.tif) DO gdal_translate -of vrt -expand rgba %%N D:\Karten\gdal\gdal2tiles\NL50\%%~nN.vrt gdalbuildvrt -allow_projection_difference ...

1

subprocess.check_output can catch anything gdal_rasterize might write to stdout or stderr. That's as good as you'll be able to do. http://docs.python.org/2/library/subprocess.html#subprocess.check_output

1

With a bit more work I've answered my own question. It required setup and testing of a new server and then application to the existing server. What I had done previously with GDAL ... enabled Ubuntu GIS "unstable" repository, installed what looked like the required update to 1.10.0: sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable sudo apt-get ...

1

Or you can use geotifcp (http://www.remotesensing.org/geotiff/geotifcp.html). To dump the data from a GeoTiff to a world file try listgeo (http://www.remotesensing.org/geotiff/listgeo.html).

1

Both your GDB file and MySQL DB do now have the spatial data -- lon-lat points -- in shape columns, but encoded as binary data. In computing, data types often have different internal and external representations. You wish to see the data presented in lon-lat, but they are represented internally in binary. You just need to use the appropriate function to ...

1

I would use a two step approach to 1) reassign raster values <3, for example, to NoData and 2) draw a polygon around the data portions of the raster dataset. For the first part, you can use Set Null (Spatial Analyst) with a statement like the following: OutRas = SetNull(InRas1 >= 3, InRas1) For the second step you can use Raster Domain (3D ...

1

This looks like a follow up on your previous post. I still advocate using 'nearblack' as I suggested there. I would do all the mosaicing, and nearblack etc in GDAL and leave converting to ECW as a final step. If you have a lot of rasters to do then you could script all of this in one script by importing both ArcPy and GDAL modules. Have a function for ...

1

the image look like there is not much dark areas inside. So you could start by computing the intensity of your RGB image (band 1 + band 2 + band 3) e.g. using raster calculator in ArcGIS, then you set a threshold (e.g. Con (intensity < 10 , 1, 0) using raster calculator). You can also do more complex testing such as directly Con(("band1" < 5) and ...

1

A first try of mine: +proj=omerc +lat_0=-21.73061814 +lonc=140.4800753 +alpha=-358.2111111 +k=0.9997418930000001 +x_0=16000 +y_0=20000 +gamma=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs using the parameters you gave, but not yet perfect. For the Affine implementation, see also: PROJ.4 custom projection that is Transverse Mercator with ...

1

You can instantiate an output driver (e.g PNG/JPEG) and use that driver to create a copy of the wms on file. I don't code in C#, but the following python snippet should help, the methods will be the same/similar. ds=gdal.Open("<GDAL_WMS>etc...</GDAL_WMS>", gdal.GA_ReadOnly) outfile='C:/Temp/test.png' pngdrv=gdal.GetDriverByName('PNG') ...

1

PosgresSQL 9.3.2.3 shouldn't be using GDAL, since it is just a database without no dependencies on it. PostGIS 2.1.1.1 doesn't require the latest GDAL (recommends at least 1.9), but if you installed it using OSGEO4W, then it is probably linked to the latest GDAL (1.10). QGIS 2.0.1 is probably linked to GDAL 1.10. QGIS 1.7 is probably linked to GDAL 1.8 or ...

1

If you add both layers to QGIS and switch On-the-fly-reprojection on, they should align. Maps in degrees (EPSG:4326) always look squeezed the more to the poles you come. So nothing wrong here. The first coordinate pair are metres, measured from 27°E +500km and the equator. The second ones are degrees. They are not rectangular anymore due to the ...

1

For Linux is something like #!/bin/bash for FILE in arcgis_folders_here/* do echo "Processing $FILE file..." gdal_translate -a_srs "+proj=tmerc +lat_0=39.668258 +lon_0=-8.133108 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs" -of GTiff$FILE \$FILE.tif done just replace "+proj=tmerc +lat_0=39.668258 +lon_0=-8.133108 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 ...

1

I had this same problem and needed to batch rectify TIFFs georeferenced in ArcGIS 10.0. First I tried to produce GeoTiffs with QGIS 2.2.0 (Raster, Conversion, Translate) which uses GDAL. First problem: GDAL will create a GeoTIFF, but it does not have any projection. This appears to be a bug. The ESRI world file in ArcGIS has an extension of *.tfwx (if ...

Only top voted, non community-wiki answers of a minimum length are eligible