New answers tagged

0

The image/raster coordinates of a point x,y (world coordinates) are gt = (263104.72544800001, 10.002079999999999, 0.0, 155223.647811, 0.0, -10.002079999999999) x,y = (263220.5,155110.6) rasterx = (x - gt[0]) / gt[1] rastery = (y - gt[3]) / gt[5] print rasterx,rastery 11.5750475901 11.3024301945 This allows to extract the value of the pixel at coordinate ...


2

I found the solution (maybe this will help someone), I have to write the Spatial reference into the CRS variable, then netcdf can recognize all projection info which is defined. crso.spatial_ref = str(srs) str(srs) gives me: crs#spatial_ref=PROJCS["Azimuthal_Equidistant", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS ...


1

Some things to think about. If your example is any indicator, you are going to have over 5 million records when you are finished. Neither shapefiles or personal geodatabases will handle that very well, so you should use either a file geodatabase or SDE. Processing 48 million records won't be quick, especially when it involves a spatial selection. Perhaps ...


2

I believe your second solution is way slower that the first one, since Clip has an inner Selection By Location plus other geometrical operations. So I suggest some improvements to your first algorithm: As per @artwork21 comment, merge the shapefiles if possible. e.g. they share the same schema. This procedure removes an extra loop over the shapefiles, ...


6

I think the speed problem you are having is looping around each feature and using geoprocessing tools inside the loop. They are not designed for that. They expect to process the whole dataset at once. So restructure your script to avoid the loop and it should complete "within the time to drink a cup of coffee". The data is in my opinion so small that it ...


0

As noted in the comments, 102020 is not an EPSG code, but an ESRI one. So the safe solution is to use gdalwarp -ot Float32 -s_srs EPSG:4326 -t_srs "+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -r near -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 etopo1.tif etopo1_warped_cli.tif In Desktop GIS like QGIS a ...


1

Although quite old, I stumbled upon this page with the same question- How do I fill the extents of a raster to a defined tile grid? The problem I was having is that all grid squares at the edge of my raster that were partially filled by the raster extents were being missed by a gdal_translate loop. Seeing as it took me a while to figure out the issues I was ...


4

I don't understand your problem import osgeo.ogr as ogr import osgeo.osr as osr driver = ogr.GetDriverByName("ESRI Shapefile") srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # create the data source data_source = driver.CreateDataSource("test_ogr.shp") # create the layer layer = data_source.CreateLayer("test_ogr", srs, ogr.wkbPoint) # Add the field ...


1

No Data values are not stored as "No Data" inside raster files but like any other value. The difference is that the dataset then has a flag which tells QGIS which value to ignore. You need to find out what value represents No Data in your raster and enter that value into the No Data field in the proximity analysis settings. You can find the value either ...


1

Use band.XSize instead of width e.g. fmt = "<" + ("h" * band.XSize) In the ReadRaster use row and column numbers, not the real coordinates e.g. for y in range(band.YSize): scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType) Check the data type with gdalinfo if it is really 2byte integer.


0

The projection looks like it matches this custom entry in http://spatialreference.org: PROJCS["Stereographic_North_Pole", GEOGCS["GCS_Coordinate System imported from GRIB file", DATUM["D_unknown", SPHEROID["Sphere",6371229,0]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], ...


-1

In case you need gdal+python on windows x64, you should compile it manually.


1

Put the proj4 string between quotes. Now the string gets split at parsing and GDAL thinks that the +ellps is a name of the source file.


2

Using the EPSG works: gdalwarp -t_srs EPSG:4269 RipBuf100_09.tif RipB uf100WGS_09.tif Creating output file that is 55104P x 24192L. Processing input file RipBuf100_09.tif. Using internal nodata values (e.g. -2.14748e+009) for image RipBuf100_09.tif. Copying nodata values from source RipBuf100_09.tif to destination ...


1

shot in the dark, but maybe use the EPSG instead of the actual CRS label


1

The correct answer was given yesterday on IRC in the GDAL channel. The above *.nc file is a netCDF version 4 encapsulated into a HDF5 format. Under Linux the libraries support the version 4 of netCDF, so the file opens with this (specific) driver and the medatadata (as min/max values) are correctly displayed. Under Windows (at least with GDAL installed ...


2

I'd check the following: first of all, I'd like to be sure that source data are OK, by comparing them with official/cadastral data I'd try to use the qgis openlayers plugin to check if the transformation to 3857/900913 works properly for both 23031 and 25831 Try to reproject using QGIS, just saving the shapefile with the target srid. If using QGIS works ...


1

I don't know why the order of drivers is different between Windows and Linux: Windows: gdalinfo --formats | findstr "netCDF HDF" HDF4 (ros): Hierarchical Data Format Release 4 HDF4Image (rw+): HDF4 Dataset HDF5 (ros): Hierarchical Data Format Release 5 HDF5Image (ro): HDF5 Dataset netCDF (rw+s): Network Common Data Format Linux: gdalinfo --formats|grep ...


0

Would go for GDALwarp as well. Be sure to be consistent with "postings" and "cell" interpretations of rasters. http://www.remotesensing.org/geotiff/spec/geotiff2.5.html


2

Your values aren't in 0,255 since they are UInt16. You can try rescaling to 0,255 (GDAL works it out by default from input min/max and output default 0,255): gdal_translate -b 1 -b 2 -b 3 -mask "none" "input.tif" "output.tif" -scale Note you can add params if the defaults aren't sensible: -scale [src_min src_max [dst_min dst_max]] ...


2

Try clipping the polygons before using them (also, please try to provide complete code including library calls in the future): library(ggmap) library(rgdal) library(rgeos) library(ggplot2) URL <- ...


0

I've marked Reilly answer as the answer as that is what I am using, but thought I'd go into a bit more detail as to how to actually go about doing this. So to actually get the filename from Python to C/C++ you need to pull in the python dataset(this doesn't have much overhead as it just takes a reference). then use PyObject* fileString = ...


2

Whenever I need to process large amounts of data I run gdal from the command line. Check http://www.gdal.org for info, examples and a couple tutorials to get you comfortable using it from the command line.


1

So you just want to mosaic all the tiles for a given day? That's a perfect Job for GDALs VRTs. gdalbuildvrt mosaic_049.vrt 049*.hdf or from Python for all days import subprocess import glob for day in range(0, 365): day = str('%0.3d' % day) cmd = ["gdalbuildvrt", "mosaic_"+day+".vrt", glob.glob(day+"*.hdf"] subprocess.call(cmd) edit: ...


2

Solved using Programmatic raster-vector calculation I had some troubles using directly gdal.RasterizeLayer() on layers like propose here in the cookbook but it seems that using "MEM" data source and finally writing it on the disk is maybe better. My solution (which is one of many possible solutions) def rasterizer(shapePath, rasterPath, attribute, ...


1

Are you sure that the driver is implemented? I have no experience with node-gdal myself, but I don't see it in the list of bundled drivers. There is also an open issue concerning PostgreSQL support. Otherwise you might want to add a driver to the gdal.open() method, but as I said, I couldn't find it in the bundled drivers. I looked for the PG: string which ...


1

Tidiest approach that I have seen is to use python to trigger the pgsql2shp command line utility which is part of the postgis package. Not very pythonic, but very simple and robust. e.g. See How can I get a shapefile from a postgis query? Once you have a shapefile it is straight forward to convert to most common formats using ogr2ogr, also over the ...


0

If you're only reading the dataset, why pass the data structure to C++? Why not pass the filename and let C++ open it shared? That said, I think the C# swig wrapper has a method to get the unmanaged (i.e. C++) handle of the dataset. It is a member of the dataset object itself. See if the Python wrapper has the same thing. Although this is more dangerous ...


0

I would reproject the files with gdalwarp. I've done the same for files in EPSG:3763 that I want to convert to EPSG:3857. I compared the results using QGIS and Geoserver and the generated images were fine. Since a small rotation is applied to the images, you might get some black lines on the border (but these lines can be made transparent afterwards). ...


0

The problem is that with easy_install GDAL, the script try to first install the last version of GDAL (2.0.1) with many errors (answer of Kersten) while what you want are the Python bindings of GDAL 1.11 Therefore, you have to look for and download the Python bindings of GDAL 1.11 and after CFLAGS=`/path of/gdal-config --cflags` LDFLAGS=`/path of ...


0

If you want to compile GDAL from source you will need to fulfill all its dependencies. Unless you explicitly need to compile it from source I would advise to use a pre-built binary. On OsX the most popular (and probably easiest) way to install GDAL is via the kyngchaos.com binaries. Alternatively if you use (Ana)Conda you can also install GDAL binaries ...


2

You could use another option like: os.system("gdal_translate -of XYZ " + source_file + " " + out_file)


5

You are trying to execute a command line utility from within Python. To do this you can use subprocess, which takes the command line arguments as a list of strings. import subprocess cmd = ["gdal_translate", "-of", "XYZ", "dataset.tif", "D:\Delft3D\Teste\GDAL\clc_raster.tif.xyz"] subprocess.call(cmd) edit: Additionally I'd like to point out that a ...


0

I ended up building a function that uses a reference image with the correct projection and extent to transform the XYZ ascii file into a GeoTiff. My ascii file had no header, so gdal_translate didn't work. Here's the function: ```def ascii_to_tiff(infile, outfile, refIm): """ Transform an XYZ ascii file without a header to a projected GeoTiff ...


1

You have Georeferenced images that you want to Rectify. There doesn't seem to be a command line Tool that mimics the one on the Georeferencing Toolbar, Rectify. I found a script, by Rob, that reads the control Points from the .tif.aux.xml file and passes this info to the Warp command, outputing a warped image. Please see the script at the end of this post, ...


0

The sd flag will retrieve a subdataset that's inside the file (if present). It's explained in gdalinfo man page. The number after the sd flag will be the subdataset number. But be carefull. For MOD13Q1, NDVI is subdataset 1. Subdataset 2 is EVI. And I believe the same order applies to MYD13A2 You can see the subdataset numbers by issuing gdalinfo on the ...


1

You are trying to adress a subdataset inside a HDF container directly. There are two ways you can do that with gdalinfo: Put the complete name of the subdataset in parantheses gdalinfo " HDF4_EOS:EOS_GRID:MYD13A2.A2015297.h16v05.005.2015314081208.hdf:MODIS_Grid_16DAY_1km_VI:1 km 16 days NDVI" Use the subdataset option gdalinfo -sd 1 ...


5

Linux (ubuntu as you call it) uses a shell (probably bash in this case) which has a different syntax to windows so you want something like: for i in se70*.tif do gdalinfo $i done


4

Instead of doing the reclassification as a double for loop described by dmh126, do it using np.where: # reclassification lista[np.where( lista < 200 )] = 1 lista[np.where((200 < lista) & (lista < 400)) ] = 2 lista[np.where((400 < lista) & (lista < 600)) ] = 3 lista[np.where((600 < lista) & (lista < 800)) ] = 4 ...


2

you can use gdalbuildvrt instead of gdalmerge to create a virtual raster template (xml file) that will behave like a merged raster. then you can run gdal_translate with the -projwin option for your set of 256*256 tiles. For most of my applications, I would create a single output per face then use gdalbuildvrt for virtual tiles. Note that if you create one ...


2

You should be able to just use the sql argument in ogr2ogr. For instance, with the following polygon Shapefile with two features and two attributes: $>ogrinfo -so -al polygon.shp INFO: Open of `polygon.shp' using driver `ESRI Shapefile' successful. Layer name: polygon Geometry: Polygon Feature Count: 2 Extent: (-1.206294, -0.828671) - (0.727273, ...


2

You should exchange ymin and ymax in -te (target extent). -te 102.3375307079206 10.350077240783799 107.6323189079206 14.6874037407838


3

The spatial reference of the tif is not defined. You will need to define it before proceeding with the clip. C:\Users\frogman\Downloads>gdalinfo Hillshade_Wrap.tif Driver: GTiff/GeoTIFF Files: Hillshade_Wrap.tif Size is 9600, 15120 **Coordinate System is ** Origin = (1773760.000000000000000,5897760.000000000000000) Pixel Size = ...


0

Using the core GDAL utilities, you could easily convert your ASCII file to GTiff using gdal_translate as follows: gdal_translate -of 'GTiff' input.asc output.tiff


0

The answer to the question is there's no need to convert the data to a bitmap format at run time. You can pass the data obtained by the RasterIO function to your opencv functions. This question is for when I worked in a company as an apprentice. The time when I didn't know anything about GDAL. So I've asked the question generally and it did not receive much ...


0

tilt pointed out in the comments above that the rastbandarg version of ST_MapAlgebra's callback form accepts more than two raster bands. This approach solves the problem.


1

ArcGIS supports a number of methods to georeference an image. IMHO its a matter of taste, what methods are associated with georeferencing, and what with "undoing" projection or distortion. The simplest of them, the affine or 1st Order Polynomial supports translation, rotation, and scaling. This one is parameterized by the numbers saved in a "world file", ...


5

You have now a better way to do. Since RFC 59.1 : GDAL/OGR utilities as a library, you can use gdalwarp from Python directly without using any call to the command line utility but using really the function from Python. This solution is a bit "on the edge" as you need at the moment to use the latest GDAL version (version 2.1, in fact the master/trunk ...


3

To convert an ASCII file with longitude, latitude and data value you may use a function like this: from osgeo import gdal def csv2tif(source, target): cvs = gdal.Open(source) if cvs is None: print 'ERROR: Unable to open %s' % source return geotiff = gdal.GetDriverByName("GTiff") if geotiff is None: print 'ERROR: ...


2

If QGIS asks you for a CRS, it does not write that information into the file. You have to save the data into a new file, then QGIS writes the CRS information into it. Alternatively, use gdal_translate -a_srs to assign the CRS to the output file. If you need a different projection, use gdalwarp -s_srs -t_srs to do both things in one step.



Top 50 recent answers are included