Re: [gdal-dev] Question on building multi band composite and going back to RGB GeoTiff

[Thomas Knudsen via gdal-dev]

Way back in a different career stage, I published this method for
displaying colour-infrared as fairly convincing pseudo-true colour:
https://www.researchgate.net/publication/242078137_A_simplified_method_for_generation_of_pseudo_natural_colours_from_colour_infrared_aerial_photos

At least one company provided material based on this algorithm, as a
commercial service around 2005.

/Thomas

Den man. 22. apr. 2024 kl. 20.13 skrev Javier Jimenez Shaw via gdal-dev <
gdal-dev@lists.osgeo.org>:

> In addition to all the previous, what about gamma correction?
> Maybe your RGB images are gamma corrected, while the other bands are not.
> In that case you will be comparing apples and oranges (in addition to the
> scaling problems described before).
> That also means that using the multispectral images (even if you had a
> Blue one) without gamma correction will have a strange color for us humans
> using a monitor.
>
> What I did in the past was to "fake" blue as a combination of red and
> green (yes, it is not perfect, but better than nothing), and apply a gamma
> of 2 to all the bands.
>
> Other "nice" option is to display them as CIR
> https://www.mngeo.state.mn.us/chouse/airphoto/cir.html
>
> Cheers
> Javier
>
> On Mon, 22 Apr 2024 at 17:49, Andrew C Aitchison via gdal-dev <
> gdal-dev@lists.osgeo.org> wrote:
>
>> On Mon, 22 Apr 2024, Raley, Nathan via gdal-dev wrote:
>>
>> > Hmm, good catch.  Looking at the stats for the red band:
>> >
>> > Band 1 Block=128x128 Type=UInt16, ColorInterp=Gray
>> >  Min=130.000 Max=36265.000
>> >  Minimum=130.000, Maximum=36265.000, Mean=10415.962, StdDev=3502.933
>> >  NoData Value=0
>> >  Metadata:
>> >STATISTICS_MAXIMUM=36265
>> >STATISTICS_MEAN=10415.961859513
>> >STATISTICS_MINIMUM=130
>> >STATISTICS_STDDEV=3502.9325925887
>> >STATISTICS_VALID_PERCENT=49
>> >
>> > So, is there no way to actually do this via command line, or am I going
>> to have to write a python lib to actually read through these and determine
>> what they should be scaled to?
>>
>> There is gdal_transform -scale (and -scale_2 etc.)
>>
>> What are you going to do with the resulting file ?
>> QGIS allows you to make these sorts of adjustments at view-time.
>>
>> Remember that the max and min values present in an image may not represent
>> the full range and the next-door image may have a larger or smaller range
>> (or just a shifted range). Imagine a row of three tiles; the middle on on
>> the edge of a glacier, the first entirely on the glacier and the last
>> entirely off the glacier. You may wish to apply the same scaling to all
>> three.
>>
>> > From: Daniel Evans 
>> > Sent: Monday, April 22, 2024 9:29 AM
>> > To: Raley, Nathan 
>> > Cc: 'gdal-dev@lists.osgeo.org' (gdal-dev@lists.osgeo.org) <
>> gdal-dev@lists.osgeo.org>
>> > Subject: [External] Re: [gdal-dev] Question on building multi band
>> composite and going back to RGB GeoTiff
>> >
>> > Hi Nathan, My initial suspicion might just be that the scaling the data
>> provider did to go from the raw data to a human-eye-friendly RGB composite
>> isn't the conversion you're assuming. I know that with the data I regularly
>> work with,
>> >
>> > Hi Nathan,
>> >
>> > My initial suspicion might just be that the scaling the data provider
>> did to go from the raw data to a human-eye-friendly RGB composite isn't the
>> conversion you're assuming.
>> >
>> > I know that with the data I regularly work with, it may be provided as
>> Uint16, but the data range doesn't extend all the way to 65535.
>> >
>> > If you compare the values in the separate R and G images to the RGB
>> composite, do they appear to match the conversion you're assuming, or is
>> there a different scaling (and possibly offset)?
>> >
>> > Cheers,
>> > Daniel
>> >
>> > On Mon, 22 Apr 2024, 15:20 Raley, Nathan via gdal-dev, <
>> gdal-dev@lists.osgeo.org> wrote:
>> > I currently have a RGB geotiff composite image that has a Byte
>> datatype.  I also have individual band images for R, G, Redge, and NIR that
>> are UInt16 datatypes.  Since I’m missing the Blue band from the individual
>> bands, I was attempting to extract the blue band from the RGB composite
>> image, scale it up to the UInt16 datatype, and build a composite VRT with
>> R, G, B, NIR, RE bands included in it.  I was then attempting to extract
>> the RGB bands from the multispec VRT in order to see if the process was
>> working as intended, but I’m getting an extremely blue image.
>> >
>> > Can anyone shed some light as to what I may be doing wrong here?
>> >
>> > I started by building a VRT for each band:
>> > gdal_translate source_RGB.tif b.vrt -ot UInt16 -of VRT -b 3 -scale 0
>> 255 0 65535
>> > gdalbuildvrt -b 1 r.vrt source_R.tif
>> > gdalbuildvrt -b 1 g.vrt source_G.tif
>> > gdalbuildvrt -b 1 nir.vrt source_NIR.tif
>> > gdalbuildvrt -b 1 re.vrt source_RE.tif
>> >
>> > I then merged the VRTs:
>> > gdalbuildvrt -separate multispec.vrt r.vrt g.vrt b.vrt nir.vrt re.vrt
>> >
>> > I now

Re: [gdal-dev] [EXTERNAL] Re: Lack of srs importFromESRI integer function

[Thomas Knudsen via gdal-dev]

In the DK registry we used descriptive names, rather than numbers to refer
to the transformations taking legacy systems to ETRS89. Essentially these
transformations amounted to CRS definitions. Is it safe to assume a
"crs_code" is strictly numeric?

Den ons. 10. jan. 2024 kl. 20.46 skrev Meyer, Jesse R. (GSFC-618.0)[SCIENCE
SYSTEMS AND APPLICATIONS INC] via gdal-dev :

> To be clear, I wasn’t advocating for less functionality.  I wanted more!
> We know what our constraints are, and we’ve had significant issues in the
> past with slightly different projection related strings producing
> surprisingly different outputs.  Keying off a single value simplifies a
> tremendous amount for us, with a slight loss in generality and increase in
> error modes.
>
>
>
> For now, the logic resembles:
>
> Crs_string = “EPSG:{crs_code}”
>
> If crs_code > 32767:
>
> Crs_string = “ESRI:{crs_code}”
>
>
>
> SetFromUserInput(crs_string)
>
>
>
> The lack of such a function implied that there were overloaded values
> between the authorities (and perhaps, there were more than 2).
>
>
>
> Jesse
>
>
>
> *From: *Javier Jimenez Shaw 
> *Date: *Wednesday, January 10, 2024 at 1:55 PM
> *To: *Even Rouault 
> *Cc: *"Meyer, Jesse R. (GSFC-618.0)[SCIENCE SYSTEMS AND APPLICATIONS
> INC]" , "gdal-dev@lists.osgeo.org" <
> gdal-dev@lists.osgeo.org>
> *Subject: *Re: [gdal-dev] [EXTERNAL] Re: Lack of srs importFromESRI
> integer function
>
>
>
> *CAUTION:* This email originated from outside of NASA.  Please take care
> when clicking links or opening attachments.  Use the "Report Message"
> button to report suspicious messages to the NASA SOC.
>
>
>
> This is the discussion I have with my colleagues every time (but with
> EPSG). The wanted to use an integer, and I ask them again and again to use
> a better notation like "EPSG:{code}". The method that Even mentions,
> SetFromUserInput, works with "EPSG:code", with "ESRI:code", with a WKT,
> with a PROJJSON, etc. It makes the management of the CRS much easier. When
> you need something more complicated than an ESRI code, if your API is only
> ready for an integer value... you have a lot of things to refactor (months
> of work, I can tell you).
>
>
>
> Technically ESRI and EPSG there are overlapping. As you can see in this
> link there are some non-deprecated ESRI codes below 33000:
>
>
> https://crs-explorer.proj.org/?ignoreWorld=false&allowDeprecated=false&authorities=ESRI&activeTypes=PROJECTED_CRS,GEOGRAPHIC_2D_CRS,GEOGRAPHIC_3D_CRS,GEOCENTRIC_CRS,GEODETIC_CRS,VERTICAL_CRS,COMPOUND_CRS&map=osm
>
>
>
> And many more if you "Allow deprecated".
>
>
>
> On Wed, 10 Jan 2024 at 19:01, Even Rouault via gdal-dev <
> gdal-dev@lists.osgeo.org> wrote:
>
>
>
> Le 10/01/2024 à 18:52, Meyer, Jesse R. (GSFC-618.0)[SCIENCE SYSTEMS AND
> APPLICATIONS INC] a écrit :
>
> Thanks Even, we’ll try that.
>
>
>
> Are there any known integer values that ESRI and EPSG conflate?  If not,
> it would be convenient to simply pass in an integer to a single API
> function and be done with it.
>
> I believe not. Codes in the [1,32767] range should lead to equivalent
> definitions under both authorities. But when they are really the same, to
> save space and confusion to users, they are not imported under the ESRI
> authority and are only available under the EPSG ones. So codes in the
> [1,32767] under the ESRI authorities are basically for CRS names that are
> slightly different spelled by ESRI (but otherwise with a definition
> compatible of the EPSG one). There are also very ancient deprecated CRS
> that were imported from EPSG, but were deprecated by EPSG so long ago that
> they are not in the EPSG database anymore (like
> https://spatialreference.org/ref/esri/31491/).
>
> But there's no only EPSG and ESRI in life :-) The IAU authority for
> example uses numeric code in the ranges of what could be used by EPSG or
> ESRI.
>
> So it is really the tuple (authority_name,code) that conveys a primary key.
>
>
>
> Jesse
>
>
>
> *From: *Even Rouault 
> 
> *Date: *Wednesday, January 10, 2024 at 12:48 PM
> *To: *"Meyer, Jesse R. (GSFC-618.0)[SCIENCE SYSTEMS AND APPLICATIONS
> INC]"  ,
> "gdal-dev@lists.osgeo.org" 
>  
> *Subject: *[EXTERNAL] Re: [gdal-dev] Lack of srs importFromESRI integer
> function
>
>
>
> *CAUTION:* This email originated from outside of NASA.  Please take care
> when clicking links or opening attachments.  Use the "Report Message"
> button to report suspicious messages to the NASA SOC.
>
>
>
> Jesse,
>
> You can use SetFromUserInput("ESRI:")  . ImportFromEPSG(code) more or
> less does SetFromUserInput("EPSG:{code}")
>
> Even
>
> Le 10/01/2024 à 18:44, Meyer, Jesse R. (GSFC-618.0)[SCIENCE SYSTEMS AND
> APPLICATIONS INC] via gdal-dev a écrit :
>
> Hi,
>
>
>
> Our team has moved away from providing geographic / projected CRS codes as
> strings and prefer using integers where possible.  However, in the C++ API,
> I only see an integer based importFromEPSG function, which doesn’t appear
> to accept ESRI codes (crs n

Re: [gdal-dev] Requiring numpy for the Python bindings

[Thomas Knudsen via gdal-dev]

Den tirs. 5. dec. 2023 kl. 02.15 skrev Greg Troxel via gdal-dev <
gdal-dev@lists.osgeo.org>:
>  Building the
> rust compiler requires the *immediately preceding* compiler, which seems
> unprecented, and I don't understand how people can think that's ok

I don't think this is true - it's merely the common case:
Over at https://rustc-dev-guide.rust-lang.org/building/bootstrapping.html
the description is

> The stage0 compiler is usually the current beta rustc compiler and its
associated
> dynamic libraries, which x.py will download for you.
> (You can also configure x.py to use something else.)

which makes it hardly any different than any compiler bootstrapping process.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Cannot find proj.db

[Thomas Knudsen]

Elena,

If your environment variables are set exactly as you write:

> PROJ_LIB= “D:\MDT\V9\bin\gdal\gdal_latest\share”
> PROJ_DATA= “D:\MDT\V9\bin\gdal\gdal_latest\share”

you should leave out the apostrophes and the extra space after the equals
sign, i.e.

set PROJ_LIB=D:\MDT\V9\bin\gdal\gdal_latest\share
set PROJ_DATA=D:\MDT\V9\bin\gdal\gdal_latest\share

The PROJ documentation on the subject is over at
https://proj.org/usage/environmentvars.html

/Thomas

Den tors. 30. mar. 2023 kl. 14.01 skrev Elena Ruiz :

> Hello,
>
> I am using GDAL 3.6.2 on Windows 11 and am having problems with the
> following error:
>
>
>
> ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db
>
>
>
> I have set as environment variable PROJ_LIB=
> “D:\MDT\V9\bin\gdal\gdal_latest\share”  and PROJ_DATA=
> “D:\MDT\V9\bin\gdal\gdal_latest\share”, as user variable and in the system
> variable, but it still doesn't work and the file proj.db exits in
> D:\MDT\V9\bin\gdal\gdal_latest\share  .
>
>
>
> D:\MDT\V9\bin\gdal\gdal_latest\x64>gdalwarp -s_srs EPSG:4326 -t_srs
> EPSG:25831 -r near -of GTiff
> C:\Users\elena\Downloads\D37531385_0101_DTM.asc
> C:/Users/elena/Downloads/D37531385_0101_DTM_25831_qgis_.tif
>
> ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db
>
> ERROR 1: Translating source or target SRS failed:
>
> EPSG:4326
>
>
>
> The funny thing is that I have the same version of GDAL installed on my
> computer using Anaconda3, to do tests, and from the anaconda directory it
> works for me, and if I look in the environment variables I don't see any
> path that refers to where Anaconda has the Proj.db. Is there another way to
> set the PROJ_LIB path besides adding an environment variable? Is there any
> way to know where the PROJ_LIB is looking for my version of GDAL?
>
>
>
> C:\Users\elena\Anaconda3\envs\mygdal\Library\bin>gdalwarp -s_srs EPSG:4326
> -t_srs EPSG:25831 -r near -of GTiff
> C:\Users\elena\Downloads\D37531385_0101_DTM.asc
> C:/Users/elena/Downloads/D37531385_0101_DTM_25831_qgis_.tif
>
> Creating output file that is 1003P x 1201L.
>
> Processing C:\Users\elena\Downloads\D37531385_0101_DTM.asc [1/1] : 0Using
> internal nodata values (e.g. -) for image
> C:\Users\elena\Downloads\D37531385_0101_DTM.asc.
>
> Copying nodata values from source
> C:\Users\elena\Downloads\D37531385_0101_DTM.asc to destination
> C:/Users/elena/Downloads/D37531385_0101_DTM_25831_qgis_.tif.
>
> ...10...20...30...40...50...60...70...80...90...100 - done.
>
>
>
> Thank you and regards
> --
>
> *Elena Ruiz *
> Sofware Development & Technical Support
> Tel. +34 952 43 97 71
> er...@aplitop.com
> Sumatra, 9 - 29190 Málaga (Spain)
>
>
> *www.aplitop.com 
>    
> 
> *
>
> In accordance with the provisions of the European Regulation of Data
> Protection 2016/679 (Reglamento Europeo de Protección de Datos 2016/679),
> we inform you that the data and the information you provide us through this
> medium will be used by APLITOP, S.L., with C.I.F. B-92543396 and with
> address at C / Sumatra, 9, Malaga, 29190, in order to answer your questions
> and inform you about our products. The data provided will be kept as long
> as it does not request its cessation and will not be transferred to third
> parties except in cases where there is a legal obligation. You have the
> right to access your personal data, correct inaccurate data or request its
> deletion when the data is no longer necessary for the purposes that were
> collected, as well as any rights recognized in the RGPD 2016/679
> 
>
>
> 
> 
>
> 
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] why padfExtent/s in gdal_alg.h ?

[Thomas Knudsen]

>  Ah ok, mostly I was just confused why it's not an error for the compiler

The compiler only needs to know the type of the argument - in a header
file, the identifier name following
the type is for human information only

Den ons. 8. feb. 2023 kl. 00.49 skrev Michael Sumner :

>
>
> On Wed, Feb 8, 2023 at 10:44 AM Even Rouault 
> wrote:
>
>> Michael,
>>
>> Feel free to submit a fix as there's obviously no justification for the
>> discrepancy other than just being accidental.
>>
>
> Ah ok, mostly I was just confused why it's not an error for the compiler
> ... but I'm getting used to surprises learning C++ :)
>
> PR incoming, thanks!
>
> Cheers, Mike
>
>
>
>> There are likely hundreds of similar occurrences in the code base.  I
>> believe there's a compiler warning (or at least cppcheck) to detect this,
>> but it is not activated.
>>
>> Even
>> Le 08/02/2023 à 00:35, Michael Sumner a écrit :
>>
>> This is not causing any problems afaict, but I happened to notice the
>> plural for 'padfExtents' in the gdal_alg.h signature for
>> GDALSuggestedWarpOutput2 compared to gdaltransformer.cpp
>>
>>
>> https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdal_alg.h#L264
>>
>>
>> https://github.com/OSGeo/gdal/blob/579693c2975cda414af75239e09b8b7b952029ca/alg/gdaltransformer.cpp#L412
>>
>>
>> The discrepancy is visible in the doc here:
>>
>>
>> https://gdal.org/api/gdal_alg.html#_CPPv424GDALSuggestedWarpOutput212GDALDatasetH19GDALTransformerFuncPvPdPiPiPdi
>>
>> Cheers, Mike
>>
>>
>>
>> --
>> Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> Hobart, Australia
>> e-mail: mdsum...@gmail.com
>>
>> ___
>> gdal-dev mailing 
>> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>> -- http://www.spatialys.com
>> My software is free, but my time generally not.
>>
>>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Create tile index by time?

[Thomas Knudsen]

Sounds like a STAC related task to me - https://stacspec.org/en

Den fre. 14. okt. 2022 kl. 09.49 skrev Stefan Gofferje <
li...@home.gofferje.net>:

> Morning all!
>
> I have a little hobby project where I pull Sentinel imagery of my town's
> area and create various composites. I run that every morning, so I have
> a bunch of raster images which represent the same geographic area at
> different times. I recently also started playing with mapserver and I
> would like to set up a mapserver to serve those images. I know mapserver
> basics, have my first mapserver up and running.
>
> What I'm looking for now is how to create a tile index for a bunch of
> images which includes the times of those images, so I can server them
> from mapserver as a layer with a time dimension.
>
> Any pointers, tutorials, "RTFMs" (with links) would be very much
> appreciated!
>
> -Stefan
>
>
> --
>   (o_   Stefan Gofferje| SCLT, MCP, CCSA
>   //\   Reg'd Linux User #247167   | VCP #2263
>   V_/_  https://www.gofferje.net   | https://www.saakeskus.fi
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdaldem aspect is grid-oriented

[Thomas Knudsen]

The proj_factors(...) function, which a.o. returns the meridian convergence
in the PJ_FACTORS typedef, may be useful in this work.

Den fre. 25. jun. 2021 kl. 12.26 skrev Even Rouault <
even.roua...@spatialys.com>:

> Michael,
>
> A general solution would have more appeal, although I suspect it will be a
> bit more computation intensive as you'll probably have to transform to
> invoke PROJ multiple times per grid point to compute the direction to the
> North
>
> Even
> Le 25/06/2021 à 06:29, Michael Sumner a écrit :
>
> Hello, we're using `gdaldem aspect` on polar grid products and we need to
> convert the grid-oriented directions to compass oriented ones,
> location-specific directions, where "compass north" is dependent on
> longitude.
>
> It's easy but inefficient for us to do this after running gdal dem for our
> polar grids at latitude -90, and a more general solution would require a
> general approach to "where north is" for any location in a projection.
>
> I'll be able to modify GDAL to do this as an option to gdaldem for a
> simple polar grid (e.g. EPSG:3031), and that will make our use of it much
> more efficient. Is there any interest in a more general approach?
>
> Best, Mike
>
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsum...@gmail.com
>
> ___
> gdal-dev mailing 
> listgdal-dev@lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Geotransform for non-North-up images?

[Thomas Knudsen]

Perhaps you can find some relevant information on the Wikpedia page
describing the ESRI world file format:
https://en.wikipedia.org/wiki/World_file - that usually where I go when I
need to re-learn the basics of affine transformations

/thomas

Den man. 21. dec. 2020 kl. 12.34 skrev G Seiffert :

> Hi all from GDLA.
>
> Thanks for the opportunity to as a question. It's regarding the Geotransform
> Tutorial (https://gdal.org/tutorials/geotransforms_tut.html). Tried to
> get info in the web but since this seems a tricky one, my searches failed.
>
> The tutorial only deals with the ideal case of 'North up' images, for
> which GT(2) and GT(4) are zero. However, my images are 'East up' and
> potentially 'any direction up' (underwater photomosaic surveys by ROV). The
> rotation works with '90' for GT(2, 4). But the scaling seems completely
> ignored. Any hint? If it would be easy, I assume your tutorial would give
> an example for how to deal with "non-N-up' images, but ...
>
> My pics are 1920x1080, with pixel resolution of 0.0015 (yes, 0.15cm per
> pixel, we're flying just 3m above the bottom).
>
> In case I rotate the images prior to geotransform (in Photoshop),
> geotransform works perfect, with GT(2, 4) = 0. Scaling spot on. I can live
> with that for our last survey but I'm also looking for a solution in case
> our survey heading cannot be 0, 90, 180, or 270, but has to be something
> like 35° (due to bottom currents etc.).
>
> Best regards, any hint qappreciated,
> Gerhard
>
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Ellipsoidal length of a line

[Thomas Knudsen]

Awesome!

Den lør. 15. jun. 2019 kl. 23.04 skrev Even Rouault <
even.roua...@spatialys.com>:

> $ cat test.csv
> id,WKT
> 1,"LINESTRING(2 49,3 50,4 49)"
>
> $ ogr2ogr test.db test.csv -f sqlite -dsco spatialite=yes -a_srs EPSG:4326
>
> $ ogrinfo test.db -sql "select st_length(geometry, 1) from test" -al -q
>
> Layer name: SELECT
> OGRFeature(SELECT):0
>   st_length(geometry, 1) (Real) = 265450.955822012
>
>
> Cf http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html#p6
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Ellipsoidal length of a line

[Thomas Knudsen]

In the PROJ package, the geod utility is included - you could do

echo lat0 lon0  lat1 lon1 | geod -I +ellps=GRS80

Where the "-I" indicates that you want to calculate the "inverse geodetic
problem", i.e. you know where you are, and where you want to go, but you
need to know how far to go and in which direction.

e.g.

$ echo 55 12 56 12 |geod -I +ellps=GRS80
> 0d -180d 111332.699

Where the distance (in m) you look for is the third number on the output
line. The two first are the forward and return azimuths respectively.
Otherwise, you could use the "proj_lp_dist()" entry of the PROJ API - cf.
https://proj.org/development/reference/functions.html

(and probably you should continue this thread on the PROJ mailing list if
you need additional assistance)

/Thomas

Den lør. 15. jun. 2019 kl. 06.43 skrev Nicolas Cadieux <
nicolas.cadi...@archeotec.ca>:

> Thanks,
> Could work but I think this will be too slow. I wonder how QGIS does it? I
> guess they use code from Proj.4.  If anyone has an other idea, shoot!
> Cheers
> Nicolas
>
> Le 14 juin 2019 à 13:20, Patrick Young 
> a écrit :
>
> Not exactly what you want, but you can do this with PostGIS by casting
> your geometry to the geography type:
>
> https://postgis.net/workshops/postgis-intro/geography.html
>
> On Fri, Jun 14, 2019 at 10:26 AM Atle Frenvik Sveen 
> wrote:
>
>> Hi!
>>
>> Not sure if would want to use gdal for this task*, but take a look at
>> this blog post:
>>
>> https://janakiev.com/blog/gps-points-distance-python/
>>
>>
>> *or if it's doable, i guess not, since the scope of gdal is
>> reading/writing geospatial formats
>>
>> -a
>>
>> -
>>   Atle Frenvik Sveen
>>   a...@frenviksveen.net
>>   45278689
>>   atlefren.net
>>
>> On Fri, Jun 14, 2019, at 17:32, Nicolas Cadieux wrote:
>> > Hi,
>> >
>> > I am trying to get the length of a line in python. (Not just the
>> > straight length between the first and last nodes).  Using geopandas,
>> > (therefore the Shapely lib) I am getting the euclidien distance even
>> > though the dataframe holdings the line geometries has a CRS (WGS84,
>> > zone UTM 18 S).  Obviously, the WGS84 Ellipsoid is not taken into
>> > account.
>> >
>> > Can I do this with gdal/ogr?
>> >
>> > Thanks for the help
>> > Nicolas
>> > ___
>> > gdal-dev mailing list
>> > gdal-dev@lists.osgeo.org
>> > https://lists.osgeo.org/mailman/listinfo/gdal-dev
>> ___
>> gdal-dev mailing list
>> gdal-dev@lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Oracle spatial with gdal

[Thomas Knudsen]

try user:password instead of user/password.

(Also, you *may* need a port number, server:port/service_name)


Den fre. 1. feb. 2019 kl. 13.56 skrev san619 :

> Thanks for the reply jmckenna.I got the driver working now.But when im
> trying
> to connect to the oracle im getting error "ERROR 1: ORA-12170: TNS:Connect
> timeout occurred"
> What could be the reason.I have installed oracle instant client and set the
> oracle path also.But im unable to connect why? Im using the query "ogrinfo
> OCI:/@/" in the gdal sdk shell
> but im getting error. i tried other way also
> "ogrinfo -ro OCI:/@/" but still the
> same error what could be the reason.
>
>
>
> --
> Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Perforce issue in coordinate transformations with Gdal 2.3.2

[Thomas Knudsen]

Prepare/finalize functions were added to pj_inv & pj_fwd  in PROJ PR #731 (
https://github.com/OSGeo/proj.4/pull/731), in order to resolve ticket #700 (
https://github.com/OSGeo/proj.4/issues/700).
These add some overhead as well but IIRC, more like +15% (much depending on
the actual projection) than +100%.

Den man. 28. jan. 2019 kl. 17.38 skrev Even Rouault :

> Stefan,
>
> I bet you also upgraded your PROJ version in the process ? I can't think
> of a reason in GDAL itself for a performance change. But on the PROJ side,
> yes, there might be some explanation
> The main one I can see is that in PROJ 4.9.3, the Extended Transverse
> Mercator is used by default for UTM, which has probably a higher
> computation time. If you define the GDAL configuration option /
> environmenet variable OSR_USE_ETMERC to NO, then the older Transverse
> Mercator method will be used.
> PROJ 5.0 has also undergone major changes that might have increased a bit
> the computation time, but I wouldn't expect them to cause a 2x slowdown.
>
> Even
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Kriging Interpolation

[Thomas Knudsen]

Hi Ian,

If you find IDW on N points slow, you will find Kriging on the same number
of points unbearably slow, since for Kriging you will have to solve a
system of N linear equations (essentially inverting an NxN matrix) to
compute the Kriging weights. To set up the matrix you will even need to
compute the same distance elements you need for IDW. So Kriging will never
become faster than IDW.

Hence, the only realistic way to use Kriging on large data sets is to base
each estimate on local subsets - typically by working in a TIN model, just
like in the linear interpolation case. And if your data set is dense, the
difference between Kriging and linear interpolation will be miniscule.
Kriging will, however, provide a variance estimate for the prediction,
which may be realistic assuming you provide a realistic covariance function
for the phenomenon under investigation.

So to make any difference wrt the linear interpolation, you will need to
include a larger number of data, e.g. by including neighbouring triangles
of the TIN, or using another data structure to subset your observations.

Some years ago I wrote an informal note about this subject. The title is
“Data Assimilation for Updates of Digital Terrain Models” and it is
available over at http://sdfe.dk/media/2916626/kms_technical_report_14.pdf

/Thomas




2018-06-25 5:19 GMT+02:00 Ian Reese :

> Hi Gdal,
>
> I'd like to know if it is possible to make Kriging Interpolation part of
> gdal_grid.  We work with a good deal of scattered point elevation data and
> neither IDW or Linear interpolations(smooth or not smooth)  produce the
> desired visual results.  Nor are they fast enough for the density of data
> we commonly work with. Linear interpolation is a reliable and fast method,
> however the triangulated results are difficult to work with visually.
>
> For some visual comparisons and reference, here is a link to a folder of
> sample points(shp) for a test area, two examples of kriged data for that
> test region, and a BASH script to reproduce IDW and Linear interpolations
> of the same region.
>
> https://drive.google.com/open?id=1HKbK4wMmVN6JZSvhxBQDub1WXBLE710H
>
> Here is an example of the Kriging method used by ArcGIS:
> http://desktop.arcgis.com/en/arcmap/latest/extensions/
> geostatistical-analyst/kriging-in-geostatistical-analyst.htm
>
> Few graphs explaining best interpolation methods to use with given data.
> Gridded vs. scattered.
> https://www.neonscience.org/spatial-interpolation-basics
>
> If it is possible, we are curious to know what type of funding would be
> needed.
>
> Cheers,
>
> Ian
>
> --
>
> This message contains information, which may be in confidence and may be
> subject to legal privilege. If you are not the intended recipient, you must
> not peruse, use, disseminate, distribute or copy this message. If you have
> received this message in error, please notify us immediately (Phone 0800
> 665 463 or i...@linz.govt.nz) and destroy the original message. LINZ
> accepts no responsibility for changes to this email, or for any
> attachments, after its transmission from LINZ. Thank You.
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev