Tag Info

New answers tagged

0

GDALbuildVRT http://www.gdal.org/gdalbuildvrt.html has a switch -vrtnodata value [value...]: this allows you to set multiple input nodata values, provided both values are OK to be set as NoData (i.e. -32768 and 255 for all images) otherwise you will have to translate somehow so they all have the same nodata value. GDAL_Translate with -a_nodata value will set ...


0

If you do not need to have a physical mosaic file, but rather the results of the mosaic, you could use a Virtual Raster (Raster - Misc - Build Virtual Raster). With this tool you can generate a mosaic'd image from rasters, and identify no data values.


0

You are running into a memory error. According to this mailing list posts from the developers: http://lists.maptools.org/pipermail/fwtools/2006-October/000531.html http://lists.osgeo.org/pipermail/gdal-dev/2010-January/023310.html you should better use gdalwarp for such purposes.


0

Modelbuilder! Add iterator: FOR Set it to start at 10, iterate from 0-350 or so, in iterations of 30. Output of iterator will include Value variable. Select by attribute, as so: FID = %Value%+1 or FID = %Value%+2 or FID = %Value%+3 or etc. For the FIDS: First iteration: 0+1, 0+2, 0+3 Second iteration: 30+1, 30+2, 30+3 Third interation: 90+1, 90+2, ...


0

Python IDE which is installed with ArcGiS can't use GDAL modules. You have few options: Install QGis and use python shell in QGis, Install new Python IDE and again install GDAL libraries, Install development enviroment like Eclipse, there you can set python interpreter and connect libraries installed somwhere else. I think that first option will be ...


0

I advise you to uninstall Qgis and after to do a fresh installation from here, all dependencies (gdal,proj4, geos, etc) will be automatically satisfied. I use Linux Mint 13 "Maya", so for me Ubuntu 12.04 (Precise) works fine. For you I think it will be Ubuntu 13.10 (Saucy).


1

The simple approach is to loop through each polygon, filter the remaining polygons by the spatial extent of the polygon (plus a little buffer), and then run your adjacency test. If your polygons are going to be of the "Simple Feature" variety, i.e. you aren't using a topological data model, you'll need to consider how you define adjacency in terms of ...


4

When you start with a Python module, there are several solutions to find the available functions. One of them is dir: geom = feat.GetGeometryRef() print dir(geom) ['AddGeometry', 'AddGeometryDirectly', 'AddPoint', 'AddPoint_2D', 'Area', 'AssignSpatialReference', 'Boundary', 'Buffer', 'Centroid', 'Clone', 'CloseRings', 'Contains', 'ConvexHull', 'Crosses', ...


0

Answering the comment about the world file, the world file doesn't explicitly say which projection is being used by the tiff, though it does encapsulate the information. So yes you can directly use the world file information, or you can run gdalinfo and it should provide the projection information for you, it's just a convenience. For more information on ...


2

It doesn't look like your code properly saves/closes the dataset. To do this, add this to the end: dst_ds = None # save, close Also, although it looks like you want to use -999 for NODATA, this needs to be set to the resulting band. If you want to learn more about raster processing with Python, check out rasterio.


0

Without having gone through your code just a suggestion: NDVI is an index and results are between 0 and 1. you will probably work with 8-bit Tiff which stores values between 0 and 255 so if you multiply your results by 100 it should work ndvi = np.where ( check, (nir - red ) / ( nir + red ) * 100, -999 )


0

I made a version of WorldFileTool (which I'll call 0.3.8 because the project seems dead) that supports CLI input and can be included as part of a GDAL script: 0.3.8 WorldFileTool jar package download link (personal build) Modified source code tarball (based on subversion repository) Command line argument: java -jar WorldFileTool.jar ...


1

Hmm as I recall there should be in an sfcgal output in your postgis_full_version() output. Did you compile postgis with sfcgal support? It's not enough to just have sfcgal installed. I see you are right the instructions in docs don't tell you how to compile with sfcgal support. I'll amend that. What you need to add is in your postgis configure ...


3

If you want to do this "right" (taking into account the fact that latitude and longitude are angular units, and using an ellipsoid as a model of the earth's shape), you can try using the geographiclib library, which is a Python version of Charles Karney's Algorithms for geodesics. See also the Wikipedia page Geodesics on an ellipsoid for a look into some of ...


4

It is easiest with shapely: from shapely.geometry import box extents = [(-180.0, -90.0, 180.0, 83.624), (-124.731, 24.956, -66.97, 49.372), (-122.42, -37.818, 151.207, 52.516)] for i in extents: a = box(i[0],i[1],i[2],i[3]) print i, a.area (-180.0, -90.0, 180.0, 83.623999999999995) 62504.64 (-124.73099999999999, 24.956, -66.969999999999999, ...


6

Main problem is getting the area of the extent. I wrote a quick ogr function to do this def extentArea(extent): #Unpack extent tuple to coordinates minX, minY, maxX, maxY = extent #unpack the tuple #Create empty geometry and add vertices geom = ogr.Geometry(type = ogr.wkbLinearRing) geom.AddPoint_2D(minX,minY) ...


0

GDAL works with GCP, but that requires a further step: Apply the GCP Transform the picture using the GCP to the source CRS Reproject from the source CRS to your desired CRS To apply the GCP correctly, you have to know the local coordinate systems used for native tif and your source CRS. The tif has the origin in the upper left corner, X to the right, and ...


0

Gene's answer pointed me in the right direction. The simple answer seems to be that the data in PostGIS is in Web Mercator (EPSG 3786) and so if I want lat/long it needs to be re-projected (or de-projected...), with an SRID like EPSG 4326. call (["ogr2ogr", "-f", "GeoJSON", record[1] + ".json", 'PG:dbname=\'gis\'', "-sql", 'SELECT ...


2

OpenLayers uses the EPSG:3857 coordinate system, in meters, and not the WGS84 system, in degrees, look at OpenStreetMap Wiki: EPSG:3857 But why use subprocess and ogr2ogr? 1) you can use directly the PostGIS ST_AsGeoJSON function: import psycopg2 conn = psycopg2.connect("dbname='osm' host='localhost' user='me'") cur = conn.cursor() # srid of the layer ...


2

For the Mercator projection, the extent can not reach North and South pole for mathematical reasons. The standard Google and Openstreetmap mercator projection is limited to 85.011° North and South to get a square map. See http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#X_and_Y for explanation. Using EPSG:3857, the extent of a map is ...


2

Individual bands can be accessed by calling GetRasterBand(4) from your datasource. You could then write your band as array into a newly created copy. For instance like this: driver = gdal.GetDriverByName("....") tDs = driver.Create(output, cols, rows, 1, gdal.GDT_Float32) ds_in = gdal.Open('in.tif') array = ds_in.GetRasterBand(4).ReadAsArray() # get ...


0

I am slowly working on a gdal app that calculates viewshed with an eye on performance. It is still in development, and no where near done, but it may be useful. It'd actually be nice to have some input, as I just did it for fun. You can find it on my github account(see below). The nearest neighbor sampling is pretty much a joke, but it's fast. Linear is ...


1

Go to the Python GDAL/OGR Cookbook 1.0 documentation and you'll have the answers to all your questions: from osgeo import ogr driver = ogr.GetDriverByName('ESRI Shapefile') shape = driver.Open('my.shp') layer= shape.GetLayer() # the crs crs = layer.GetSpatialRef() and you can also create a projection file if the shapefile does not ...


1

ogrinfo layer.shp layer -so will give you general information about the shapefile.


0

In ModelBuilder, you use Iterate Feature Classes to loop through the shapefiles, then use Collect Values to send the iterator values to Merge. It will look like this:


7

It would be easier with arcpy. for i in range(10000): pol_list = [] for j in range(30): pol_list.append("a" + str(i*30 + j + 1) + ".shp") arcpy.Merge_management(pol_list, "b" + str(i+1) + ".shp") EDIT: for a feature class inside a geodatabase, you don't need the + ".shp" anymore, and you can define the workspace using : ...


0

The problem is that I was not creating a field to store the raster band. After digging through the gdal_polygonize.py file, I realized this is not automatically done when calling gdal.Polygonize, which instead uses the function found here. Here is the extra step needed to create a field and write a band to the field: newField = ogr.FieldDefn('MYFLD', ...


3

gdaldem will use the colors of the lookup table at their exact position and it will blend those colours by interpolating between the input values. so, in your example, the value of 40 will be coloured exactly in 140 0 3 and the value of 35 coloured in 204 9 9, then the "missing" values are interpolated (you can see on the answer that you mentioned that you ...


1

You could use the OTB applications : first rescale your data in 8 bit otbcli_Rescale -in input_image.tif -out intermediate_image.tif uchar -outmin 0 -outmax 255 then apply a colormap otbcli_ColorMapping -in intermediate_image.tif -method custom -method.custom.lut yourLUT.txt -out output_image.jpg your LUT should go from 0 0 255 0 to 255 ...


1

You could use the R "raster", "rasterVis" and "ggplot2" packages to automate this quite nicely. You have considerable control of very simple (just the raster) or customized plots using the R low-level plotting engine or higher level plotting like ggplot2. You can also easily call other plotting devices to output other formats including: tiff, bmp, png or ...


4

Evil Genius seems to have a good idea, but actually the color-relief file can be set using percentages as well as min and max. So you could write the file using something like: 0% green 100% red So the command would be something like: gdaldem color-relief inputfile.img colorfile.txt output.jpg -of "JPEG2000"


3

You need to have admin privileges on the box. For Windows 7: Click on the start menu. Right click on "Computer" Click "Properties" Click "Advanced system settings on the left menu bar You may have to type in a password here. The second box is titled "system variables" on of them is "PATH" Select path, click edit. Go to the end of the PATH variable ...


3

A couple of options: Convert the PostGIS layer to a shapefile in QGIS (save-as), then use the Vector|Conversion|Rasterize tool; Use the gdal_rasterize command directly. For the second option: gdal_rasterize -a VAL -ts [x] [y] PG:'host=localhost dbname=DB user=USER' -sql "SELECT the_geom, VAL FROM table" out.tif Where: VAL = the value to assign ...


1

You could try the Rasterize tool in QGIS (Raster - Conversion)


0

You can use GDAL to convert each band to an array and get the mean of the array with Numpy. import gdal, numpy ds = gdal.Open("yourRaster.tif") bandCount = 1 while bandCount <= ds.RasterCount: band = ds.GetRasterBand(bandCount) bandArray = band.ReadAsArray() print numpy.mean(bandArray) bandCount += 1


1

255 is the default NoData value in QGIS. I am not sure, what exactly the problem is with the way you tried it, but you could use the GDAL Python bindings to do what you want. For instance the following script converts your shp to a polygon based on the attribute NAME_2_NUM. Import the libraries import ogr, gdal, osr Open your shapefile source_ds = ...


0

I think your first option, memory allowing, is to read all the images into an Ndarray, tuning your storage type as @sgillies says - especially if you're just making graphs. Feed your grapher slices of the array, not new dictionaries. At the other end of the spectrum, if you're really memory constrained, is to read individual pixel values for each x and y ...


0

There's no ultimate answer to this because there's no getting away from memory vs speed tradeoffs. One way to get 4x "more" memory would be to convert your NDVI arrays from float to signed Int16. You'd still have 4 decimal places precision, and a lot of the float precision was false anyway since the source imagery would have been (I presume) 8 or 16 bit to ...


1

You don't need GDAL's Python module to do this, you can use the gdal_translate program to subset images: http://www.gdal.org/gdal_translate.html. See the -srcwin and -projwin options.


0

I have no problem with the speed of the calculation with the function rasterize: library(raster) ## Set up a raster #Extent ext <- extent(munisCrop.shape) #Resolution xy <- abs(apply(as.matrix(bbox(ext)), 1, diff)) r <- raster(ext, ncol=xy[1]/0.1, nrow=xy[2]/0.1) ## Rasterize the shapefile #you need to define the the value(s) to be transferred ...


1

I'd recommend removing the old GDAL completely. As an alternative, just remove the binaries that cause problems (probably from /usr/local/bin with something like sudo rm /usr/local/bin/gdal_translate. If you need to use both versions, you can probably set LD_LIBRARY_PATH to select particular library variants and use explicit paths to the binary you want. ...


1

I'm not sure why you're accessing img_stack[y, 0, z] instead of img_stack[x, y, z] But let's make it easier - if you want the maximum values in a numpy array along your z axis, you should be able to use: final_values = numpy.amax(img_stack, axis=2) see numpy.amax


1

Be wary of stating, out the gate, that something is producing "wrong" results. The phenomena of slope percent approaching infinity as slope degrees approach 90 is well known. You could just truncate slopes > 100 == 100. ESRI actually provides a very nice description on how slope is calculated. You could try the calculation in degrees to make sure that ...


1

I had to do something rather like this for my masters thesis, but with much fewer observer and target points. I'm not aware of a reasonable way to create a complete raster of "visible area," at least not one that wouldn't take a long time. Repeatedly running Viewshed, once for each centroid of the raster's cells, would certainly work... but as you've, ...


0

Files with non-ASCII names are not supported as you can see here: https://github.com/anitagraser/TimeManager/issues/40


0

You can turn off the error messages with gdal.PushErrorHandler('CPLQuietErrorHandler'). If you want to log the errors to a file, you can use a try: except: with gdal.UseExceptions() as per @nickves' answer, or you can pass a custom error function with gdal.PushErrorHandler(mycustomfunction).


0

You have a problem with your PATHS: with Homebrew, GDAL is installed in /usr/local/Cellar with symbolic links of the programs (ogr2ogr, etc.) in /usr/local/bin with the Frameworks of KyngChaos, the versions of GDAL are installed in /Library/Frameworks/GDAL.framework/Versions/1.x/GDAL with the programs in ...


1

You may want to enjoy the new temporal GIS framework in GRASS GIS 7: GRASS as Temporal GIS presentation PDF A temporal GIS for field based environmental modeling (article) Manual pages: http://grass.osgeo.org/grass70/manuals/temporalintro.html An initial release of GRASS GIS 7 has been done two days ago at the Vienna OSGeo Code sprint: ...


2

I recommend using Python or R (or a GIS software), as @Marc Pfister has suggested. However, you can do it with bash and gdal only, and heavy usage of grep. First get the Min/Max values without coordinates: Obtain the Min / Max values with gdalinfo or gdalinfo -mm like explained in your other question about Min/Max values. Use grep (and possibly some awk) ...


4

Is Python an option? Use RasterIO (a Python GDAL/ numpy bridge) to load the raster to NumPy array, then use numpy.amax() to find the maximum value, followed by numpy.where() to find the row/column indices, then calculate the lat and lon from the raster extents.



Top 50 recent answers are included