Vjylku

Most bit efficient text communication method?

The Nth Gryphon Number

Sliceness of knots

How to write capital alpha?

Karn the great creator - 'card from outside the game' in sealed

Why can't I install Tomboy in Ubuntu Mate 19.04?

Why are vacuum tubes still used in amateur radios?

A term for a woman complaining about things/begging in a cute/childish way

Misunderstanding of Sylow theory

How does the math work when buying airline miles?

How can I prevent/balance waiting and turtling as a response to cooldown mechanics

What makes a man succeed?

What would you call this weird metallic apparatus that allows you to lift people?

Do I really need to have a message in a novel to appeal to readers?

Lagrange four-squares theorem --- deterministic complexity

Did Mueller's report provide an evidentiary basis for the claim of Russian govt election interference via social media?

What are the discoveries that have been possible with the rejection of positivism?

Can a Beast Master ranger change beast companions?

Why does 14 CFR have skipped subparts in my ASA 2019 FAR/AIM book?

Co-worker has annoying ringtone

Antipodal Land Area Calculation

How were pictures turned from film to a big picture in a picture frame before digital scanning?

How long can equipment go unused before powering up runs the risk of damage?

What is the difference between a "ranged attack" and a "ranged weapon attack"?

gdal_translate with custom corner coordinates in python

.everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty{ margin-bottom:0;
}

0

I wish to convert this Python code using a gdal command line tool to using the gdal python library. I am using it to convert some PNG satellite images I am downloading to georeferenced TIFFs. I manually specify the bounding box of the images when I download them (in EPSG:3857), so I know all the info about the images, I just wish to store all the info in a single georeferenced TIFF file.

My current solution is to use gdal_translate called by subprocess but apparently this can all be done natively in python, so I would prefer to do that rather than rely on this hacked-together workaround.

from osgeo import gdal

import subprocess

import shlex    



# ll_x etc... are all corner coords I calculate before hand 

subprocess.run(shlex.split("gdal_translate -of GTiff -a_ullr {} {} {} {} -a_srs {} {} {}".format(

    ll_x, ur_y, ur_x, ll_y, 

    "EPSG:3857", 

    "input_sat_image.png", 

    "referenced_sat.tif")))

I have seen individual questions that cover pieces of this, but I am struggling to put them all together.

python gdal gdal-translate

share|improve this question

edited Aug 6 '18 at 20:31

PolyGeo♦

54k1782246

asked Aug 6 '18 at 20:17

dirtshelldirtshell

112

bumped to the homepage by Community♦ 3 mins ago

This question has answers that may be good or bad; the system has marked it active so that they can be reviewed.

add a comment | 

0

I wish to convert this Python code using a gdal command line tool to using the gdal python library. I am using it to convert some PNG satellite images I am downloading to georeferenced TIFFs. I manually specify the bounding box of the images when I download them (in EPSG:3857), so I know all the info about the images, I just wish to store all the info in a single georeferenced TIFF file.

My current solution is to use gdal_translate called by subprocess but apparently this can all be done natively in python, so I would prefer to do that rather than rely on this hacked-together workaround.

from osgeo import gdal

import subprocess

import shlex    



# ll_x etc... are all corner coords I calculate before hand 

subprocess.run(shlex.split("gdal_translate -of GTiff -a_ullr {} {} {} {} -a_srs {} {} {}".format(

    ll_x, ur_y, ur_x, ll_y, 

    "EPSG:3857", 

    "input_sat_image.png", 

    "referenced_sat.tif")))

I have seen individual questions that cover pieces of this, but I am struggling to put them all together.

python gdal gdal-translate

share|improve this question

edited Aug 6 '18 at 20:31

PolyGeo♦

54k1782246

asked Aug 6 '18 at 20:17

dirtshelldirtshell

112

bumped to the homepage by Community♦ 3 mins ago

This question has answers that may be good or bad; the system has marked it active so that they can be reviewed.

add a comment | 

I wish to convert this Python code using a gdal command line tool to using the gdal python library. I am using it to convert some PNG satellite images I am downloading to georeferenced TIFFs. I manually specify the bounding box of the images when I download them (in EPSG:3857), so I know all the info about the images, I just wish to store all the info in a single georeferenced TIFF file.

My current solution is to use gdal_translate called by subprocess but apparently this can all be done natively in python, so I would prefer to do that rather than rely on this hacked-together workaround.

from osgeo import gdal

import subprocess

import shlex    



# ll_x etc... are all corner coords I calculate before hand 

subprocess.run(shlex.split("gdal_translate -of GTiff -a_ullr {} {} {} {} -a_srs {} {} {}".format(

    ll_x, ur_y, ur_x, ll_y, 

    "EPSG:3857", 

    "input_sat_image.png", 

    "referenced_sat.tif")))

I have seen individual questions that cover pieces of this, but I am struggling to put them all together.

python gdal gdal-translate

share|improve this question

edited Aug 6 '18 at 20:31

PolyGeo♦

54k1782246

asked Aug 6 '18 at 20:17

dirtshelldirtshell

112

from osgeo import gdal

import subprocess

import shlex    



# ll_x etc... are all corner coords I calculate before hand 

subprocess.run(shlex.split("gdal_translate -of GTiff -a_ullr {} {} {} {} -a_srs {} {} {}".format(

    ll_x, ur_y, ur_x, ll_y, 

    "EPSG:3857", 

    "input_sat_image.png", 

    "referenced_sat.tif")))

share|improve this question

edited Aug 6 '18 at 20:31

PolyGeo♦

54k1782246

asked Aug 6 '18 at 20:17

dirtshelldirtshell

112

share|improve this question

edited Aug 6 '18 at 20:31

PolyGeo♦

54k1782246

asked Aug 6 '18 at 20:17

dirtshelldirtshell

112

edited Aug 6 '18 at 20:31

PolyGeo♦

54k1782246

edited Aug 6 '18 at 20:31

PolyGeo♦

54k1782246

asked Aug 6 '18 at 20:17

dirtshelldirtshell

112

asked Aug 6 '18 at 20:17

dirtshelldirtshell

112

1 Answer
1

active

oldest

votes

0

I might be wrong, but maybe you can try this way :

from osgeo import gdal 

import osr



inputImage = YOURNAME 

outputImage = YOURNAME



dataset = gdal.Open(inputImage) 

I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)



outdataset = gdal.GetDriverByName('GTiff') 

output_SRS = osr.SpatialReference() 

output_SRS.ImportFromEPSG(3857) 

outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0],gdal_Float32) 

for nb_band in range(I.shape[0]):

    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])



gcp_list = [] 

gcp_list.append(gdal.GCP(lon,lat,alt,col,row)) #Repeat this step for as many GCP as you want (there is a limit but I don't know the number of GCPs)



outdataset.SetProjection(srs.ExportToWkt()) 

wkt = outdataset.GetProjection() 

outdataset.SetGCPs(gcp_list,wkt)



outdataset = None

It worked for me, but it was in a specific case (and I did not need a very precise result...). I think you might have a problem with the sizes, but I actually don't know how to solve it.

share|improve this answer

answered Aug 7 '18 at 9:13

ArthurWillArthurWill

64

add a comment | 

Your Answer

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "79"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: false,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});


}
});

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name

Email

Required, but never shown

StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fgis.stackexchange.com%2fquestions%2f291912%2fgdal-translate-with-custom-corner-coordinates-in-python%23new-answer', 'question_page');
}
);

Post as a guest

Name

Email

Required, but never shown

0

I might be wrong, but maybe you can try this way :

from osgeo import gdal 

import osr



inputImage = YOURNAME 

outputImage = YOURNAME



dataset = gdal.Open(inputImage) 

I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)



outdataset = gdal.GetDriverByName('GTiff') 

output_SRS = osr.SpatialReference() 

output_SRS.ImportFromEPSG(3857) 

outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0],gdal_Float32) 

for nb_band in range(I.shape[0]):

    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])



gcp_list = [] 

gcp_list.append(gdal.GCP(lon,lat,alt,col,row)) #Repeat this step for as many GCP as you want (there is a limit but I don't know the number of GCPs)



outdataset.SetProjection(srs.ExportToWkt()) 

wkt = outdataset.GetProjection() 

outdataset.SetGCPs(gcp_list,wkt)



outdataset = None

It worked for me, but it was in a specific case (and I did not need a very precise result...). I think you might have a problem with the sizes, but I actually don't know how to solve it.

share|improve this answer

answered Aug 7 '18 at 9:13

ArthurWillArthurWill

64

add a comment | 

0

I might be wrong, but maybe you can try this way :

from osgeo import gdal 

import osr



inputImage = YOURNAME 

outputImage = YOURNAME



dataset = gdal.Open(inputImage) 

I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)



outdataset = gdal.GetDriverByName('GTiff') 

output_SRS = osr.SpatialReference() 

output_SRS.ImportFromEPSG(3857) 

outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0],gdal_Float32) 

for nb_band in range(I.shape[0]):

    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])



gcp_list = [] 

gcp_list.append(gdal.GCP(lon,lat,alt,col,row)) #Repeat this step for as many GCP as you want (there is a limit but I don't know the number of GCPs)



outdataset.SetProjection(srs.ExportToWkt()) 

wkt = outdataset.GetProjection() 

outdataset.SetGCPs(gcp_list,wkt)



outdataset = None

It worked for me, but it was in a specific case (and I did not need a very precise result...). I think you might have a problem with the sizes, but I actually don't know how to solve it.

share|improve this answer

answered Aug 7 '18 at 9:13

ArthurWillArthurWill

64

add a comment | 

I might be wrong, but maybe you can try this way :

from osgeo import gdal 

import osr



inputImage = YOURNAME 

outputImage = YOURNAME



dataset = gdal.Open(inputImage) 

I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)



outdataset = gdal.GetDriverByName('GTiff') 

output_SRS = osr.SpatialReference() 

output_SRS.ImportFromEPSG(3857) 

outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0],gdal_Float32) 

for nb_band in range(I.shape[0]):

    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])



gcp_list = [] 

gcp_list.append(gdal.GCP(lon,lat,alt,col,row)) #Repeat this step for as many GCP as you want (there is a limit but I don't know the number of GCPs)



outdataset.SetProjection(srs.ExportToWkt()) 

wkt = outdataset.GetProjection() 

outdataset.SetGCPs(gcp_list,wkt)



outdataset = None

It worked for me, but it was in a specific case (and I did not need a very precise result...). I think you might have a problem with the sizes, but I actually don't know how to solve it.

share|improve this answer

answered Aug 7 '18 at 9:13

ArthurWillArthurWill

64

from osgeo import gdal 

import osr



inputImage = YOURNAME 

outputImage = YOURNAME



dataset = gdal.Open(inputImage) 

I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)



outdataset = gdal.GetDriverByName('GTiff') 

output_SRS = osr.SpatialReference() 

output_SRS.ImportFromEPSG(3857) 

outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0],gdal_Float32) 

for nb_band in range(I.shape[0]):

    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])



gcp_list = [] 

gcp_list.append(gdal.GCP(lon,lat,alt,col,row)) #Repeat this step for as many GCP as you want (there is a limit but I don't know the number of GCPs)



outdataset.SetProjection(srs.ExportToWkt()) 

wkt = outdataset.GetProjection() 

outdataset.SetGCPs(gcp_list,wkt)



outdataset = None

share|improve this answer

answered Aug 7 '18 at 9:13

ArthurWillArthurWill

64

answered Aug 7 '18 at 9:13

ArthurWillArthurWill

64

answered Aug 7 '18 at 9:13

ArthurWillArthurWill

64