2

I have a hard problem to find the width/height of a 2D GeoTiff file which is reprojected from a source CRS (e.g: EPSG:4326) to target CRS (e.g: EPSG:3857) without running gdalwarp. Because, the file could be large (GBs) and I don't need the projected output file, just need the width/height of it.

What could be done with gdalwarp is:

gdalwarp  -s_srs EPSG:4326 -t_srs EPSG:3857 full.tif full_3857.tif 

gdalinfo full_3857.tif 

Size is 879, 811 (width/height).

I've searched a lot and what seems to be good without doing the real projection is this gdal function: GDALSuggestedWarpOutput2 https://github.com/OSGeo/gdal/blob/106c8288e7a05f4efc1a588c5a3b2da7ec52d915/gdal/alg/gdaltransformer.cpp#L354. However, it doesn't help because my application developed in Java and it uses gdal-java (GDAL JNI) https://search.maven.org/#artifactdetails%7Corg.gdal%7Cgdal%7C1.11.1%7Cpom as library.

Unfortunately, I cannot find this GDALSuggestedWarpOutput2() from gdal-java http://gdal.org/java/overview-summary.html then cannot invoke this C++ function from my application to test.

Can anyone please give me a hint to solve this problem?

1

The above mentioned solution is in python. Since you are using java, following is the solution in java language.

    gdal.AllRegister();

    Dataset source_file = gdal.Open("mean_summer_airtemp.tif");

    String source_srs = source_file.GetProjectionRef();
    SpatialReference sSRS = new SpatialReference(source_srs);

    SpatialReference dSRS = new SpatialReference();
    //sSRS.ImportFromEPSG(4326); // Need not import, let it be the original source SRS
    dSRS.ImportFromEPSG(3857);

    String source_wkt = sSRS.ExportToWkt();
    String dest_wkt = dSRS.ExportToWkt();

    double error_threshold = 0.125;
    Dataset reproj_file = gdal.AutoCreateWarpedVRT(source_file, source_wkt, dest_wkt, gdalconstConstants.GRA_NearestNeighbour, error_threshold); 

    System.out.println("Source raster Dimension: "+source_file.GetRasterXSize()+" "+source_file.GetRasterYSize());    
    System.out.println("Reproj raster Dimension: "+reproj_file.GetRasterXSize()+" "+reproj_file.GetRasterYSize());

It will help finding the dimensions of raster without using gdalwrap.

0

To solve this problem, one can use GDAL function: AutoCreateWarpedVRT

An example of how it works in Python to get the ul_x, ul_y, lr_x, lr_y and pixel_size_x, pixel_size_y

from osgeo import gdal,osr

source_file = gdal.Open("mean_summer_airtemp.tif")

source_srs = osr.SpatialReference()
dest_srs = osr.SpatialReference()

source_srs.ImportFromEPSG(4326)
dest_srs.ImportFromEPSG(3857)

source_wkt = source_srs.ExportToWkt()
dest_wkt = dest_srs.ExportToWkt()

error_threshold = 0.125
reproj_file = gdal.AutoCreateWarpedVRT(source_file, source_wkt, dest_wkt, gdal.GRA_NearestNeighbour, error_threshold)
reproj_attributes = reproj_file.GetGeoTransform()

print reproj_attributes

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.