Tag Info

Hot answers tagged

10

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


8

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.


8

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.


7

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


7

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


6

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


6

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'): ...


6

you need to specify the type of code that you are using. Here I guess that it is EPSG C:\Program Files (x86)\FWTools2.4.7>gdalwarp -s_srs EPSG:25830 -t_srs EPSG:23030 C:\orto.e cw C:\orto_ed50.ecw


6

Nice and reproducible question. Personally, I'd expect that the reason for the difference is in the implementations of the bilinear reprojection. You can obviously look into source code for the two approaches, but I'd expect that to be a vast overkill. It appears that the R implementation introduces bigger "errors" / "changes" than the raw GDAL version ...


5

I have a vague outline of how to do this, but there's plenty that you may have to understand. NetCDF files are complex general data containers so its not always clear how to get spatial data out of them. In this case, you can get the Soil_Moisture variable and that is just a 2d matrix with no coordinate reference. If you do image(A) you should see your soil ...


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

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


4

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


4

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


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


4

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


4

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


4

Reorganise your shapefile so that one shapefile contains one feature (A,B,C in your case) only Then use a loop like for i in A B C; do gdalwarp -cutline $i.shp ... $i.tif done to create each output raster. Example of script: #!/bin/sh # "shp" - folder for shapefiles # "outputraster" - folder for output rasters cd /home/user/cliprasters/ # ...


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


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

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

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


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

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

gdalwarp has to use Create() not CreateCopy(), and AAIGrid driver doesn't have Create(). Simple solution: gdaltest 14:40:41 $ gdalwarp -of vrt rj_DEM_1978.grd junk.vrt gdaltest 14:40:41 $ gdal_translate -of AAIGrid junk.vrt junk.grd


3

Creating a mosaic from your source images takes only a few seconds if you use GDAL virtual raster as output. Read http://www.gdal.org/gdal_vrttut.html and http://www.gdal.org/gdalbuildvrt.html Often the artifacts like you have can also be avoided by taking care that the individually warped images are aligned to use a common canvas. This can be achieved by ...


3

Try adding the -et (error threshold) option with lower thresholds than the default (0.125). When I use "-et 0.01", the horizontal artifacts disappear: gdalwarp -t_srs "+proj=stere +lat_0=90 +lon_0=-45 +lat_ts=70 +ellps=WGS84 +datum=WGS84 +units=m" \ -et 0.01 \ ./UiMbqSd.withmetadata.tif ./regridded_lon0_-45.tif


3

Here is an example that does roughly what you ask for. The main parameters are the geotransform array that gdal uses to describe a raster location (position, pixel scale, and skew) and the epsg code of the projection. With that, the following code should properly georeference the raster and specify its projection. I did not test this much, but it seemed to ...


3

I don't know if it's possible to clip a raster with an other raster but you could use gdaltindex to build the shapefile with the extent of your raster. http://www.gdal.org/gdaltindex.html



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