Tag Info

New answers tagged

1

You could use the GDALAutoCreateWarpedVRT function. This function creates a virtual raster dataset that represents the projected version of your raster. If you use this function, then you can read the properties of your image in the new coordinate system from the virtual raster's geotransform parameters.


0

GDAL 1.10 added a few resampling methods which will help, see gdalwarp. In particular, the -r average method, documented as: average resampling, computes the average of all non-NODATA contributing pixels. This isn't tested, but should look something like: gdalwarp -t_srs EPSG:4326 -tr 0.5 0.66 -r average fine_one_sq_km.tif coarse_average.tif Then to ...


0

You can try reprojecting the raster data but to EPSG:102113 instead, don't use 3857. This post may help: Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL


1

Your problem comes from adressing the subdatasets wrong. If you run gdalinfo on the complete file it will display the names of the subdatasets: SUBDATASET_1_NAME=HDF5:"A2015069000500.L2_LAC_OC.nc"://geophysical_data/Kd_490 To get the information of the first subdataset you need to feed the complete name into gdalinfo gdalinfo ...


1

Assuming the string field is named text and the field with corresponding integer values is named code and field code already exists, then the following Python code will do the job. # get active layer aLayer = qgis.utils.iface.activeLayer() # get fieldindex fni_t = aLayer.fieldNameIndex('text') fni_c = aLayer.fieldNameIndex('code') # initializations ...


1

assuming that the goal is to determine an elevation with a user-supplied set of coordinates, there should be a way to do this with gdal - i use a function like this, which came almost verbatim from this answer (assumes raster is not rotated). rast = 'myIMGfile.img' mx = XCOORDINATE my = YCOORDINATE def readRastPix(rast,mx,my): src_ds=gdal.Open(rast) ...


0

Here is something that might point you in the right direction: DataSize = nRows * nCols; byte[] buff = new byte[DataSize]; // THIS IS WHERE THE DATA WILL BE ThisBand.ReadRaster(cOff, rOff, nCols, nRows , buff , nCols, nRows , 0, 0); cOff and rOff are the offsets from the upper right cell in cell coordinates. nCols and nRows need to be specified twice, ...


1

If you are on Windows you could give the precompiled wheels by Christoph Gohlke a try, he has a stable 1.11 with bindings for Python 3.4 listed. If you are on Linux your best bet is compiling from source with support for your Python version (./configure --with-python).


0

Make sure you set the full path. This error seems to occur when you only give gdal the name of the file. In addition, only use ASCII characters in file paths. The (quite common GIS) problem is that non-ASCII characters (like æ ø å) in the path name gets deleted. It looks like GDAL eats the special characters. Use ASCII paths instead. Furthermore, you can ...


0

Converts my solution to a comment. Maybe more words makes it less trivial? Regardless, all these answers may work, but the SuperOverlay format is horrible, and the quads thing is pretty limiting/crude. I reverse engineered an output from OKMap... And you could use that, but I posted a script for ArcGIS here:Exporting 3GB ArcGIS Raster to KML without losing ...


0

Apart from gdal (see other answers) and building a virtual file, you can use the OTB library. This is a C++ open source library including a large set of filters for image processing. Specifically, the otb::MultiChannelExtractROI does the trick. It is also available as an application if you want to use it directly.


3

It is nearest neighbor as written in http://lists.osgeo.org/pipermail/gdal-dev/2006-November/010619.html There is also a hint in the mail "If you want control over the resampling used, you should use gdalwarp instead."


0

If you can do from osgeo import gdal,ogr your GDAL was built with the Python bindings otherwise you wouldn't be able to import it. If you want to start a GDAL utility (gdalwarp, gdal_translate, gdal_merge) from within Python your best bet is to use Pythons subprocess module as these utilities are spawned from the command line. An example: import ...


1

It is probably possible with GDAL but not extremely easy. Everything that can be done with SQL in SpatiaLite can also be done with GDAL http://www.gdal.org/ogr_sql_sqlite.html and this query could be used as starting point: ogrinfo -dialect sqlite test.shp -sql "select attr, st_union(geometry) from test group by attr" However, the query creates one ...


0

Since the first comment turned out to be the answer here's the slightly longer version of it: For gdal to write something to disk you need to flush/close the dataset. If you add ET = None at the end of your script this should work and will write the new data to disk. You can see the same behaviour in the source code of the gdal_fillnodata.py script. To ...


0

Creating contours is an interpolation, which is why your output looks like it does. You need to reclassify using gdal_calc: http://www.gdal.org/gdal_calc.html But there is some good information here relating to reclassifying with either gdal or python: Reclassify rasters using GDAL and Python And then convert to vector using gdal_polygonize: ...


0

i may try to attempt the answers in reverse: the resolution might be somewhat specific to what you want - and as you mention, dependent on the cell size. If your coordinates were in meters - you might want to end with 10 meter cell size, so you might use cellsize = 10 ncol = int(math.ceil(xmax-xmin)) / cellsize nrow = int(math.ceil(ymax-ymin)) / cellsize ...


1

I cannot imagine that UTM is not suitable for your purposes. Are you sure you picked the correct zone for the location of your dataset? Example: http://spatialreference.org/ref/epsg/32632/ I don't know much about P.T.L. but it seems to work similarly as UTM. It is unlikely that it is included in common GIS coordinate reference system (CRS) definitions, ...


0

Use gdal.GDT_Float32 not GDT_Byte. The Byte datatype is unsigned 8 bit integer. That means it can only hold values from 0-255. If you try to include values <0 or >255 they will overflow, i.e -1 will be converted to 255, -2 to 254 and 256 to 0 and 257 to 1 etc... Your code works fine with Float32 datatype, see image below. You may be seeing "values of ...


0

Here's a command line method using gdalwarp, formulated as a Windows batch file: @echo off setlocal set index=D:\source\nts_index_250k.shp set tiles=117D 117A 116O 116P for %%a in (%tiles%) do ( @echo gdalwarp -cutline %index% ^ -cwhere "'TILE_NAME' = '%%a'" ^ -cblend 3 ^ %1 %2_%%a.tif ) endlocal goto :eof Usage crop_by_poly.bat infile.tif ...


1

If you can see your layer on a map using L.TileLayer, the problem is in the boundary polygon you pass in L.TileLayer.BoundaryCanvas. It is difficult to guess the exact problem without example sources, but you can try the following checks: Add your layer using L.TileLayer and your boundary using L.Polygon (or using L.GeoJOSN). Check that your polygon ...


3

Here's just the nuts-and-bolts from some working code: #include "gdal_priv.h" #include "gdal_alg.h" // in main()... GDALAllRegister(); // register all drivers // open your raster - format doesn't matter as all the drivers are registered GDALDataset* SourceRasterDS = (GDALDataset*) GDALOpen(Raster,GA_ReadOnly); double GeoTransform[6]; ...


1

I use a simple way to solve this using the UNIX bash. I made a script in the same path of images and ran it ("script.sh"). #!/bin/bash for i in *.tif do gdal_calc.py -A $i -B MASK.dat --calc="A/B" --NoDataValue=0 --format=ENVI --outfile=directory/$i.bin done MASK.dat is a binary image with only 10.000 values to scale in NDVI range. "i" are all the MODIS ...


2

You need to post the input raster (in UTM) information too so a true before-and-after comparison can be made. Representing data using lat/lon in a raster means using a Plate Carree-like projection and treating the decimal degree values as if they're linear measures. UTM data is often 'tilted' in comparison so data is resampled. There'll be 'no data' values ...


0

If I'm understanding correctly, couldn't you just take the extent, round it to the next round integer and then subset the original raster to the new extent? Of course you would lose some data on the edge of the raster, but that would be necessary.


1

If you use gdal_polygonize.py to extract the shape of the nodata area (or the ! nodata) area into a polygon, you can use ST_Buffer() in PostGIS with a negative parameter to shrink that area back a little and then ST_Intersection() to clip the contour lines with it.


1

It looks very strange, but warping with QGIS (which runs its embedded gdalwarp) is much faster! I was able to process 14Gb file in 70 minutes, on windows, without much resource consumption. It still was not looking like it used multiple cores, but did the work, which is great. Also, it seems the same applies to gdal_translate. Probably, they build gdal ...


0

for one particular band, here's how I create a nodata mask. import gdal import numpy as np nodata = -9999#or set it based on the raster nodata value. r = gdal.Open(rastername) raster_arr = np.array(r.GetRasterBand(band).ReadAsArray()) nodatamask = raster_arr == nodata #do your thing and end up with #a result_raster that you'd like to save ...


0

try something like this: import gdal import numpy as np area_per_pixel = 100 #example...you change to suit r = gdal.Open(rasterfile) band = 1 raster_arr = np.array(r.GetRasterBand(band).ReadAsArray()) for cover in np.unique(raster_arr): tot_num_pixels = np.sum(raster_arr == cover) area = tot_num_pixels * area_per_pixel print cover, area as ...


4

I suggest tackling this using Virtual Raster Table (.vrt) format. How the end result is to be used will determine how many steps are needed. Simplest possible case is the end product will be used by a GDAL or GDAL-aware program, create one .vrt in the desired projection and then use that in your final program: gdalwarp -t_srs wgs84 -of vrt ...


1

I believe you're after re-sampling a raster on the fly which I don't believe is possible (though I'd be happy to be proved wrong). You'll have to create a new dataset at the required sample level. In ArcGIS you can do this using the Resample tool (which you can add to a python script). This does create a new file on your disk which depending on the size of ...


0

There are two python ways to do this simply, use os.walk() to find the files and then either of create a batch file or use subprocess.Popen() to execute the command. This code contains both.. default is Batch File which can be switched by setting UsingBat = False: import os, sys, subprocess BaseInFolder = sys.argv[1] BaseOutFolder = sys.argv[2] MaskImage ...


0

Sounds like a user issue with your computational region: if you want to export the full map, be sure to set the computational region to the full map prior to export (g.region rast=your_rastermap -p). Then export. Note that r.info reports the full map min/max while those of the actual computational region can differ from that. E.g., you have a DEM of Europe ...


0

First of all, I want to thank all the support I received from you guys, to AndreJ Luke, Nathan W, Michael Miles-Stimson and user30184. Following all your suggestions finally I could successfully install ECW SDK 3.3, I follow basically what appears in the OpenStreetMap Wiki. I describe below all the steps I used, a small modification from the steps explained ...



Top 50 recent answers are included