The following approach worked pretty well.
First I build virtual raster.
gdalbuildvrt raster.vrt -srcnodata 0 -input_file_list paths.txt
paths.txt is file with following content:
Then I add a pixel function to it, as showed here https://lists.osgeo.org/pipermail/gdal-dev/2016-September/045134.html.
Pixel function is written using numpy, ...
gdal_translate needs the coordinates of the corners - I don't know where your " 729000 528000" come from - I guess the width/height? No, it should be the coordinates:
gdal_translate -a_srs "+proj=gnom +lat_0=50.008 +lon_0=14.447" -a_ullr -301500 +217500 427500 -310500 QQNDR.png located.tiff
Another change I've made is to name the output .tiff since its a ...
I downloaded the glc_shv10_DOM.Tif and get very different gdalinfo output with Origin = (-180.000000000000000,90.000000001440014) and Pixel Size = (0.008333333333400,-0.008333333333400). You appear to have stripped the georeferencing out of the tif somehow or possibly an issue with your GDAL install.
If you just want to tile internally and specify the CRS, ...
In other words, you want to create a World file from the coordinates of the 4 corners and the width and height of the image
1) You get the width and height of the image with osgeo.gdal, rasterio or any other libraries to open image files as Pillow and others.
dataset = rasterio.open('satel.tif')
rasterx = dataset.width
rastery = dataset.height
2) you ...
From the GDAL doc:
The driver does support creating new files, but the input data must be exactly formatted as a SRTM-3 or SRTM-1 cell. That is the size, and bounds must be appropriate for a cell.
The output of yourgdal_translate command is pretty clear. Your data needs to be in WGS84 (epsg 4326), be 1201x1201 (cellsize=0.00083333333) or 3601x3601 (...
Thank you very much for posting your workflow, this helped me with a similar issue I was having. In case this might be useful to somebody else, I used different python functions for my raster mosaic. In my case, the no data value for the VRT was 255 and because my data only goes from 0 to 100, I masked all the values in my VRTs greater than 100 before ...
Your existing data is compressed. You're getting an increase in file-size because the default is uncompressed and you aren't specifying a compression option (-co compress=LZW).
However, it's a longstanding issue that gdalwarp doesn't deal with compression well. The solution is to gdalwarp without compression then gdal_translate with compression.
To avoid ...
gdalinfo on your data returned below:
Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Driver: netCDF/Network Common Data Format
Size is 5000, 3000
Coordinate System is:
If your first image is in EPSG:3857, it is a square. Therefore:
gdal_translate -of Gtiff -co "tfw=yes" -a_ullr -20037508.3427892 20037508.3427892 20037508.3427892 -20037508.3427892 -a_srs "EPSG:3857" "/espg-3857.tiff" "tfw.tiff"
Because 85.0511 is the result of an approximation (20037508.3427892 is also an approximation, but it is the GDAL's ...
The division in gdal is an integer division by default. You can change this behaviour by dividing with a float. In you case, simply replace 10000 by 10000.0
gdal_calc.py --type=Float32 -A C:\z\input.tif --outfile=C:\z\output.tif --calc="A/10000.0"
Note that storing in 16bit integer is more efficient than float 32, so I should think twice before running ...
Your 1st step is wrong, as the initial raster is in EPSG:4326 it should be:
gdal_translate -a_srs EPSG:4326 -a_ullr -104.7 30.8 -103.8 29.8 MRMS_RALA_LATEST.tiff projected.tiff
Then your 2nd step will correctly reproject the raster to epsg:3857 for web mercator.
A raster is a set of cells that form a grid; each cell has a value. When you reproject a raster, you are re-drawing the grid to be aligned with a new projection. So, in the below figure, your original raster grid is shown in blue, and the reprojected grid is shown in red.
Right away you can see a problem--the grids do not align. So, for example, in the ...
Note the difference between band rendering done by QGIS and the actual raster extent and values. The values of your raster remain the same, regardless of the band rendering (Multicolor band, Singleband gray, Paletted unique/values) you select in Layer properties.
If you take your original layer (i.e. HDF4_EOS:EOS_GRID:"MCD19A2.A2018312.h24v06.006....
after some research i did the following:
gdalwarp -s_srs "+proj=eqc +R=1737400" -t_srs "+proj=ortho +lat=90 +lon_0=0 +R=1737400" input.tif output.tif
So basically what caused the error was the -s_srs definition. Since it was the moons dem I had to specify the Radius by +R=1737400. Hope this helps!
Step 1: Go to https://search.earthdata.nasa.gov/search?q=SPL2SMAP_S
Step 2: Select the tile that pops up
Step 3: Filter granules by date and region
Step 4: Select + Add to Project on right
Step 5: Select My Project button that appears
Step 6: Select More Options
Step 7: Select Customize under Select Data Access Method
Step 8: Under Reformat Output select ...
It can be achieved with a simple walk or listdir:
import os, sys # the standard imports
BaseFolder = r'c:\your\folder\with\data' # change this to match your data
for FullPath, dirs, files in os.walk(BaseFolder):
for ThisFile in files: # iterate the files
fN,fE = os.path.splitext(ThisFile) # separate file name and extension
WarpOptions has xRes and yRes arguments (output x & y resolution).
You can either create a WarpOptions object and pass it to gdal.Warp with the options=your_warp_options argument or you can just pass xRes & yRes directly to gdal.Warp as keyword arguments:
options = gdal.WarpOptions(xRes=6.5, yRes=6.5, your other args, eg. resampleAlg=gdal....
In cases like this I always find that the documentation of the GDAL utility helps me work out what the Python is expecting. So looking at the gdalwarp page leads me to suspect that dstSRS is the same as the -t_srs argument so it's looking for the target SRS (or CRS) - in your case probably the same as the input one.
With the builtin GDAL utilities you can first make the VRT with gdalbuildvrt:
gdalbuildvrt mosaic_file.vrt *.tif
And then convert to Tiff and reproject with gdalwarp:
gdalwarp -t_srs EPSG:32734 -of GTiff mosaic_file.vrt output.tif
Gdal_translate can do that with -srcwin https://gdal.org/programs/gdal_translate.html.
-srcwin Selects a subwindow from the source image for copying based on pixel/line location.
The minimal command would be
gdal_translate -srcwin 30 30 1220 900 input.tif output.tif
Instead of using gdal_translate binary you can use it also as a library in Python ...
In documentation for gdalwarp command:
there is an example where is pointed out that options in your command are precisely used for un-georeferenced images and it is also necessary a cutline as csv file. So, for georeferenced images, you only need to add or to rest 30*(Pixel Size) to xmin ymin xmax ymax values. These ...
The main problem I see is that GDAL does not recognize the CRS/Transform information from the dataset. As such, transforming the dataset is not possible as the original information cannot be detected.
So, you need to construct the CRS/Transform yourself.
Step 0: Get the CRS of the dataset
Based on the link you gave, this should be the CRS of the ...
I managed to resolve the issue. The error in my case resulted from the image not being georeferenced. So even though the vector mask was correctly referenced compared to the original image, QGIS (I believe) needs to have those points referenced in a projected coordinate system.
So, find a georeferenced raster image of your area, create GCPs in the raster ...
I did that a few days ago. I am not an expert in Python so perhaps there are other ways to achieve that. Don't forget to import all the modules you need.
First, you can convert your DEM into an array of value 1 representing the extent of your DEM:
# input (dem)
dem = dem.tif
# output (dem extent)
extent = extent.tif
# convert the dem into an array of ...
It was not clear to me at the time that each srcnodata value corresponds to a single NoData value in a given band, i.e., -srcnodata "1 2 3..." indicates that 1 is the NoData value in band 1, 2 for band 2, and so on. As one commenter suggested, a solution would be to use gdal_calc.py or a similar map algebra procedure to collapse multiple desired NoData ...
We can do one error at a time:
Warning: The target file has a 'nc' extension, which is normally used by the GMT, netCDF drivers, but the requested output driver is GTiff. Is it really what you want?
If you don't want a geotiff output, which is the default. You need to set a format so add:
It might be a vector file as well as suggested in the ...
The gdal_translate command has a -projwin argument for selecting a window in projected coordinates:
-projwin ulx uly lrx lry
This needs upper left and lower right coordinates which you get from gdalinfo:
Upper Left (-13042354.590, 3866165.781) (117d 9'41.27"W, 32d46'48.73"N)
Lower Right (-13041857.790, 3865668.981) (117d 9'25.21"W, 32d46'35.22"N)
So the ...
Gdalwarp stumbles over the nodata values in the latitude and longitude bands of the netcdf file. The related bug issue is fixed in GDAL 2.4.2: https://github.com/OSGeo/gdal/issues/1451
As a workaround, I did this:
Extract the desired band to a vrt:
gdal_translate -of VRT HDF5:"input.nc"://geophysical_data/Rrs_655 -a_nodata -32767 input.vrt
Open the file ...
Here is a bit of bash script that takes an GML/JP2 image file from Sentinel-2B Level-2 10m resolution product and corrects its header. First it extracts the existing corner-coordinates using gdalinfo (twice); then it adds 10m to the two latitude values; then it writes the new corner coordinates back into the original image header using gdal_edit.py
The image ...