Tag Info

Hot answers tagged

6

You could explicitly set the output coordinate range using the target extent option to gdalwarp (ie. "-te -180 -90 180 90") but you can also use the CENTER_LONG configuration option to force rewrapping around a new central longitude. Something like this: gdalwarp -t_srs WGS84 ~/0_360.tif 180.tif -wo SOURCE_EXTRA=1000 \ --config CENTER_LONG 0 ...


6

While I don't know why GDAL provides this overlap in functionality, be sure to set the cache for gdalwarp to make it really fast: # assuming 3G of cache here: gdalwarp --config GDAL_CACHEMAX 3000 -wm 3000 $(list_of_tiffs) merged.tiff Be sure to not define more cache than having RAM on the machine.


5

Your source coordinate system is most likely not defined in the CVS file that GDAL searches for proj4 strings. It looks like you might be able to pass the source EPSG as 3031 (from spatialreference.org) Note that it looks like your input is in a local projection. Is this clipped from a larger raster? To explicitly define the source you could just provide ...


5

If you're using GDAL 1.8.0 (which I recommend because it adds a number of useful features), you can use: gdalwarp -cutline "PG:dbname=gisdb" -csql 'select * from polytest where id=1' -crop_to_cutline -of GTiff -srcnodata -9999 -dstnodata -9999 src.tif dest.tif Note the "-crop_to_cutline" parameter. With gdalwarp, I have found it always pays to use ...


4

I would not recommend using the MODIS sinusoidal projection in analysis. It would be prudent to project your MODIS data to something a bit more tractable. You can request MODIS in a projected geographic (lat/long) coordinate system on the MODIS Golbal Subsets site . That said I have used this as my CRS for MODIS "+proj=sinu +R=6371007.181 +nadgrids=@null ...


4

Try to specify the nodata-value from your input raster and set it for the output as well. Furthermore add the option -crop_to_cutline to make exact crops. More about the options here. gdalwarp -srcnodata <in> -dstnodata <out> -crop_to_cutline -cutline INPUT.shp INPUT.tif OUTPUT.tif


4

I noticed that gdal2tiles numbers the tiles from south to north (according to the TMS specification), while Openstreetmap and others do it from north to south. For my personal use, I changed the code of gdal2tiles to get it right again. See also: ...


3

You need to use the -dstalpha option to gdalwarp e.g.: gdalwarp -cutline INPUT.shp -crop_to_cutline -dstalpha INPUT.tif OUTPUT.tif This will add an alpha band to the output tiff which masks out the area falling outside the cutline. A late answer, but hopefully it will help someone else with the same problem.


3

It looks like your version of GDAL doesn't have the definition for Spherical Mercator, AKA Web Mercator, AKA Google Mercator. You could try using EPSG 900913 which is the old code for 3857 - you only need to run gdal_translate -a_srs EPSG:900913 ... on your sources. But ideally you should get GDAL version 1.9.0


3

Regarding execution times: Lanczos is extremely slow, but other resampling algorithms produce similar (perhaps even better) results: near (11s) bilinear (17s) average (27s, a bit soft) bicubic (30s, called cubic in GDAL, text thickness sometimes odd) antialias (45s, requires pil and numpy) cubicspline (1m53s, way too soft) lanczos (11m8s, extremely ...


3

You don't necessarily need pyramids, unless you have hundreds of GB of data. For lower amounts of data a mosaic should work just fine, provided you add inner tiling and overviews in your data. Have a look here for directions: http://demo.geo-solutions.it/share/foss4g2011/gs_steroids_sgiannec_foss4g2011.pdf


3

I think its because there's no inverse Winkel-Tripel projection in Proj4. You need the inverse to project a grid forward, since the target grid is derived by inverse projection from its grid points to the source grid. If you did it the other way, by mapping each source grid to a grid square in the target projection, you wouldn't end up with a regular grid ...


3

Your horizontal units are metres (= linear units). You can not derive your z-value units (vertical units) from the coordinate system. It's like checking the coordinate system properties to try and work out if a temperature raster z-values are in centigrade, kelvin or fahrenheit. You need to go back to the metadata, if there is any. You can probably ...


3

Here are the steps how I managed to clip your Geopdf from the last question: Read the imprinted border coordinates clockwise from the map and put them into a text file like this: Nr;WKT 1;POLYGON((-74.25 41.75, -74.0 41.75, -74.0 41.5, -74.25 41.5, -74.25 41.75)) Load the text file as text delimited file, with ; as delimiter and NAD27 EPSG:4267 as ...


3

I'm agree with Nathan. You need to pythonize your whole script. So substitute your for loop with something like the following: import os, fnmatch def findRasters (path, filter): for root, dirs, files in os.walk(path): for file in fnmatch.filter(files, filter): yield file for raster in findRasters(INPUT_FOLDER, '*.tif'): ...


3

Firstly, crop the image source (coords are expressed in pixels here) with: gdal_translate -srcwin 115 18 1360 2156 2104.gif 2104_cropped.tif Then, transform the known WGS84 coordinates of the upper left and lower right corners to the "WGS 84 / World Mercator" projection (EPSG:3395): cs2cs +init=epsg:4326 +to +init=epsg:3395 -35 75 -3896182.18 ...


3

It's a well known and longstanding issue that gdalwarp doesn't deal with compression well. The solution is to gdalwarp without compression then gdal_translate with compression. To avoid two lengthy processes, gdalwarp to VRT first, it's really quick, then gdal_translate with the -co compress=lzw option. i.e. $ gdalwarp -tap -tr 30 30 -t_srs "etc..." -of ...


3

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 ...


3

I would recommend to use gdalcopyproj.py, a sample file from the GDAL repository done for this purpose as mentioned directly in the script: Duplicate the geotransform and projection metadata from one raster dataset to another, which can be useful after performing image manipulations with other software that ignores or discards georeferencing ...


2

Have you tried specifying -s_srs "EPSG:4326" instead?


2

You can use gdal2tiles with the -k option to force a KML output. That will create a SuperOverlay you can open and visualize in Google Earth. Here are some instructions you may want to reference: https://developers.google.com/kml/articles/raster


2

Basically you need to cut the raster into two parts and piece them back together with a new offset/scale. There's an example here of how to do that from [-180,180] to [0,360] with gdal_translate and the VRT driver: http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial Scan down to the "5 min tutorial" and the details are under "Virtual Files". It ...


2

gdal2tiles uses the bottom-left point of origin, not top-left, so the filename pattern is {zoom}/{x}/{tilenum-y}.png where tilenum is the number of vertical tiles for the zoom level. It's also worth noting that gdal2tiles.py does not require the source data to be pre-warped, but it does require it to be RGB or RGBA.


2

How about using gdalwarp? First copy the above well known text to 'target_srs.prf' (or whatever name you like). Then: gdalwarp -t_srs target_srs.prf inputraster outputraster Full documentation on gdalwarp is here. I think that you can safely omit the s_srs flag as gdal should read the inputraster to get the srs. Alternatively, what is your ...


2

Found this: http://osgeo-org.1560.n6.nabble.com/Unable-to-load-PROJ-4-library-libproj-dylib-td5031442.html and so added export PROJSO=/Library/Frameworks/PROJ.framework/PROJ to my .bash_profile it works! I'm still running OSX 10.7 and I haven't needed this in my profile before so why I should need it now I can't say by the whole issue with paths is ...


2

I assume you have data in EPSG:31276 MGI Balkans 6. The correct projection string is: +proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs There are some bugs in different QGIS installations which result in a projection string with only three non-zero ...


2

You have to specifiy the path to gdalwarp (in the example below I took: "C:\python32\python.exe", but you have to change it to your intallation path) and change your call function like that: call(['C:\Program Files\GDAL\gdalwarp.exe', '-t_srs ' + crs, '-dstnodata 0', '-q', '-cutline ' + mask, '-dstalpha', '-of GTIFF', input, output]). or: ...


2

For the quality of British National Grid reprojection in QGIS, see Raster incorrectly reprojected to OSGB(27700). Note that the buffer will always be in the layers CRS and units, not the CRS of the project.


2

I finally managed with this very simple and clean script, which calls GDAL from Python without importing it (as suggested, but using "Call ()" method instead of "os.system ()". I hope this could help! import os, fnmatch from subprocess import call call(["ls", "-l"]) inFolder= ...


2

Convert the file to WGS84 gdalwarp in_test.vrt out_test.vrt -t_srs "+proj=longlat +ellps=WGS84" Calcualte the bbox with GDAL in Python import gdal ds = gdal.Open('out_test.vrt') cols = ds.RasterXSize rows = ds.RasterYSize geotransform = ds.GetGeoTransform() bb1 = originX = geotransform[0] bb4 = originY = geotransform[3] pixelWidth = geotransform[1] ...



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