[gdal-dev] gdal-dev Digest, Vol 133, Issue 17
Pieter Kempeneers
kempenep at gmail.com
Tue Jun 9 08:20:19 PDT 2015
>
> Message: 1
> Date: Tue, 9 Jun 2015 14:31:48 +0000
> From: "Ronquillo, Edgar Nahum" <eronquillo at lanl.gov>
> To: "gdal-dev at lists.osgeo.org" <gdal-dev at lists.osgeo.org>
> Subject: [gdal-dev] Overlay shapefile onto Geotiff image
> Message-ID:
> <
> 2086CE7E967C1B4F9D5897BE1C3C746F3119455B at ECS-EXG-P-MB01.win.lanl.gov>
> Content-Type: text/plain; charset="iso-8859-1"
>
> Hi,
> I have been working on this issue for quite a while. I currently have a
> GeoTiff image and a shapefile. I would like to overlay the shapefile on the
> Geotiff, is this possible with Gdal? I looked into gdal_rasterize but I
> just don't know how to incorporate it with what I need. A code to start
> with would be great if possible. Or maybe point me into the right direction
> of probably other libraries that could do this.
>
Hi Edgar,
What is it exactly you want to do with the overlay? If it is just cropping
the Geotiff, you can use gdalwarp (check option -crop_to_cutline). If you
want to extract the raster pixels in the Geotiff covered by the shapefile
(all pixels, average, etc.), you can use pkextract from pktools (
http://pktools.nongnu.org/html/pkextract.html). It also has a processing
toolbox plugin for QGIS (http://plugins.qgis.org/plugins/pktools/). The
code is open source (GPLv3).
Pieter.
2015-06-09 16:39 GMT+02:00 <gdal-dev-request at lists.osgeo.org>:
> Send gdal-dev mailing list submissions to
> gdal-dev at lists.osgeo.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
> or, via email, send a message with subject or body 'help' to
> gdal-dev-request at lists.osgeo.org
>
> You can reach the person managing the list at
> gdal-dev-owner at lists.osgeo.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of gdal-dev digest..."
>
>
> Today's Topics:
>
> 1. Overlay shapefile onto Geotiff image (Ronquillo, Edgar Nahum)
> 2. Re: assigning vertical datum to image files (Newcomb, Doug)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Tue, 9 Jun 2015 14:31:48 +0000
> From: "Ronquillo, Edgar Nahum" <eronquillo at lanl.gov>
> To: "gdal-dev at lists.osgeo.org" <gdal-dev at lists.osgeo.org>
> Subject: [gdal-dev] Overlay shapefile onto Geotiff image
> Message-ID:
> <
> 2086CE7E967C1B4F9D5897BE1C3C746F3119455B at ECS-EXG-P-MB01.win.lanl.gov>
> Content-Type: text/plain; charset="iso-8859-1"
>
> Hi,
> I have been working on this issue for quite a while. I currently have a
> GeoTiff image and a shapefile. I would like to overlay the shapefile on the
> Geotiff, is this possible with Gdal? I looked into gdal_rasterize but I
> just don't know how to incorporate it with what I need. A code to start
> with would be great if possible. Or maybe point me into the right direction
> of probably other libraries that could do this.
>
> Thanks for the help
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150609/028f5e7e/attachment-0001.html
> >
>
> ------------------------------
>
> Message: 2
> Date: Tue, 9 Jun 2015 10:38:30 -0400
> From: "Newcomb, Doug" <doug_newcomb at fws.gov>
> To: Even Rouault <even.rouault at spatialys.com>
> Cc: gdal dev <gdal-dev at lists.osgeo.org>, Hermann Peifer
> <peifer at gmx.eu>
> Subject: Re: [gdal-dev] assigning vertical datum to image files
> Message-ID:
> <
> CALQGVr2pLpssKbeenSWsk4qULUs7t1WT9-G+ThQHXQ947-mTLA at mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
>
> Evan,
> Here is the output from gdalinfo wit the expanded config.
> Note the error messages on failure to load datum shift files for the
> corner coordinates.
> I'll upgrade to RC1 and check that again, but thought I would mention.
>
> gdalinfo --config GTIFF_REPORT_COMPD_CS YES D05_37_30081003_20141231.tif
>
> Driver: GTiff/GeoTIFF
> Files: D05_37_30081003_20141231.tif
> Size is 1000, 1000
> Coordinate System is:
> COMPD_CS["unknown",
> PROJCS["NAD83(2011) / North Carolina (ftUS)",
> GEOGCS["NAD83(2011)",
> DATUM["NAD83_National_Spatial_Reference_System_2011",
> SPHEROID["GRS 1980",6378137,298.257222101,
> AUTHORITY["EPSG","7019"]],
> AUTHORITY["EPSG","1116"]],
> PRIMEM["Greenwich",0,
> AUTHORITY["EPSG","8901"]],
> UNIT["degree",0.0174532925199433,
> AUTHORITY["EPSG","9122"]],
> AUTHORITY["EPSG","6318"]],
> PROJECTION["Lambert_Conformal_Conic_2SP"],
> PARAMETER["standard_parallel_1",36.16666666666666],
> PARAMETER["standard_parallel_2",34.33333333333334],
> PARAMETER["latitude_of_origin",33.75],
> PARAMETER["central_meridian",-79],
> PARAMETER["false_easting",2000000],
> PARAMETER["false_northing",0],
> UNIT["US survey foot",0.3048006096012192,
> AUTHORITY["EPSG","9003"]],
> AXIS["X",EAST],
> AXIS["Y",NORTH],
> AUTHORITY["EPSG","6543"]],
> VERT_CS["NAVD88 height",
> VERT_DATUM["North American Vertical Datum 1988",2005,
>
>
> EXTENSION["PROJ4_GRIDS","g2012a_conus.gtx,g2012a_alaska.gtx,g2012a_guam.gtx,g2012a_hawaii.gtx,g2012a_puertorico.gtx,g2012a_samoa.gtx"],
> AUTHORITY["EPSG","5103"]],
> UNIT["metre",1,
> AUTHORITY["EPSG","9001"]],
> AXIS["Up",UP],
> AUTHORITY["EPSG","5703"]]]
> Origin = (3010000.000000000000000,805000.000000000000000)
> Pixel Size = (5.000000000000000,-5.000000000000000)
> Metadata:
> AREA_OR_POINT=Area
> DataType=Generic
> Image Structure Metadata:
> COMPRESSION=DEFLATE
> INTERLEAVE=BAND
> Corner Coordinates:
> ERROR 1: failed to load datum shift file
> Upper Left ( 3010000.000, 805000.000)
> ERROR 1: failed to load datum shift file
> Lower Left ( 3010000.000, 800000.000)
> ERROR 1: failed to load datum shift file
> Upper Right ( 3015000.000, 805000.000)
> ERROR 1: failed to load datum shift file
> Lower Right ( 3015000.000, 800000.000)
> ERROR 1: failed to load datum shift file
> Center ( 3012500.000, 802500.000)
> Band 1 Block=1000x2 Type=Float32, ColorInterp=Gray
> Description = Layer_1
> Min=-3.000 Max=-3.000
> Minimum=-3.000, Maximum=-3.000, Mean=-3.000, StdDev=0.000
> NoData Value=-3.40282306073709653e+38
> Metadata:
> LAYER_TYPE=athematic
> STATISTICS_COVARIANCES=0
> STATISTICS_MAXIMUM=-3
> STATISTICS_MEAN=-3
> STATISTICS_MEDIAN=0
> STATISTICS_MINIMUM=-3
> STATISTICS_MODE=0
> STATISTICS_SKIPFACTORX=1
> STATISTICS_SKIPFACTORY=1
> STATISTICS_STDDEV=0
>
>
> Doug
>
> On Tue, Jun 9, 2015 at 10:21 AM, Newcomb, Doug <doug_newcomb at fws.gov>
> wrote:
>
> > Evan,
> > I went back and checked. The .aux.xml sidecar files with the img files
> > do not appear to have projection information.
> >
> > Doug
> >
> > On Tue, Jun 9, 2015 at 9:32 AM, Even Rouault <even.rouault at spatialys.com
> >
> > wrote:
> >
> >> Selon "Newcomb, Doug" <doug_newcomb at fws.gov>:
> >>
> >> > Thanks!
> >> >
> >> > That should work in this case to build the virtual raster. I'm still
> >> > curious about assigning the vertical datum to the rasters.
> >>
> >> Doug,
> >>
> >> I'm a bit surprised that the HFA driver reports a SRS with a VERTCS node
> >> as a
> >> child of the PROJCS node. I believe this is an incorrect WKT definition.
> >> There
> >> should rather be a COMPD_CS root note with 2 children: a PROJCS node for
> >> the
> >> horizontal CRS and a VERTCS node for the vertical CRS. Perhaps this SRS
> >> string
> >> is defined in the .aux.xml sidecar files that are apparently present ?
> >>
> >> Regarding your question, there are a few EPSG codes that correspond to
> >> aCOMPD_CS, for example 6349. Otherwise, GDAL allows to build a custom
> >> COMPD_CS
> >> with the syntax EPSG:XXXX+YYYY where XXXX is the EPSG code of the
> >> horizontal SRS
> >> and YYYY the EPSG code of the vertical SRS. In your case that should be
> >> EPSG:6543+5703 I think
> >>
> >> So you can do things like "gdal_translate in.tif out.tif -a_srs
> >> EPSG:6543+5703"
> >>
> >> If you read back out.tif, you will only get the horizontal SRS by
> >> default. This
> >> was for backward compatibility reason I think. But if you specify
> >> "--config
> >> GTIFF_REPORT_COMPD_CS YES" you will get the COMPD_CS.
> >>
> >> Even
> >>
> >> >
> >> > Doug
> >> >
> >> > On Tue, Jun 9, 2015 at 1:48 AM, Hermann Peifer <peifer at gmx.eu> wrote:
> >> >
> >> > >
> >> > > About your initial problem:
> >> > > > but gdalbuildvrt complained, informing me that
> >> > > > they were not in the same projection.
> >> > >
> >> > > What about using the below gdalbuildvrt option?
> >> > >
> >> > > Hermann
> >> > >
> >> > > -allow_projection_difference:
> >> > > (starting with GDAL 1.7.0) When this option is specified, the
> utility
> >> will
> >> > > accept to make a VRT even if the input datasets have not the same
> >> > > projection. Note: this does not mean that they will be reprojected.
> >> Their
> >> > > projection will just be ignored.
> >> > >
> >> > > Source: http://www.gdal.org/gdalbuildvrt.html
> >> > >
> >> > >
> >> > >
> >> > > On 2015-06-08 21:01, Newcomb, Doug wrote:
> >> > >
> >> > >> Hi Folks,
> >> > >> I have a directory of 783 .img format Lidar - based DEMs in the
> same
> >> > >> projection, North Carolina State Plane Feet (2011) , NAD 83 ,
> >> NVD88. I
> >> > >> was going to use gdalbuildvrt to create a single virtual image for
> >> the
> >> > >> area, but gdalbuildvrt complained, informing me that they were not
> in
> >> > >> the same projection.
> >> > >>
> >> > >> A couple of quick bash scripts/commands
> >> > >>
> >> > >> for x in *.img; do gdalinfo $x >>img_info.txt;done
> >> > >>
> >> > >> and
> >> > >> grep PROJCS img_info.txt|sort|uniq -c
> >> > >>
> >> > >> gives me the following:
> >> > >>
> >> > >> 437
> >> PROJCS["NAD_1983_2011_StatePlane_North_Carolina_FIPS_3200_Ft_US",
> >> > >> 346
> >> PROJCS["NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet_2011",
> >> > >>
> >> > >> gdalinfo gives the following for each type of file ( origin section
> >> > >> clipped out) :
> >> > >>
> >> > >> Driver: HFA/Erdas Imagine Images (.img)
> >> > >> Files: D05_37_20878102_20141231.img
> >> > >> D05_37_20878102_20141231.img.aux.xml
> >> > >> Size is 1000, 1000
> >> > >> Coordinate System is:
> >> > >> PROJCS["NAD_1983_2011_StatePlane_North_Carolina_FIPS_3200_Ft_US",
> >> > >> GEOGCS["GCS_NAD_1983_2011",
> >> > >> DATUM["NAD_1983_2011",
> >> > >> SPHEROID["GRS_1980",6378137.0,298.257222101]],
> >> > >> PRIMEM["Greenwich",0.0],
> >> > >> UNIT["Degree",0.017453292519943295]],
> >> > >> PROJECTION["Lambert_Conformal_Conic_2SP"],
> >> > >> PARAMETER["False_Easting",2000000.0],
> >> > >> PARAMETER["False_Northing",0.0],
> >> > >> PARAMETER["Central_Meridian",-79.0],
> >> > >> PARAMETER["Standard_Parallel_1",34.3333333333333],
> >> > >> PARAMETER["Standard_Parallel_2",36.1666666666667],
> >> > >> PARAMETER["Latitude_Of_Origin",33.75],
> >> > >> UNIT["Foot_US",0.30480060960121924],
> >> > >> VERTCS["NAVD_1988_Foot_US",
> >> > >> VDATUM["North_American_Vertical_Datum_1988"],
> >> > >> PARAMETER["Vertical_Shift",0.0],
> >> > >> PARAMETER["Direction",1.0],
> >> > >> UNIT["Foot_US",0.3048006096012192]]]
> >> > >>
> >> > >>
> >> > >> Driver: HFA/Erdas Imagine Images (.img)
> >> > >> Files: D05_37_20957301_20141231.img
> >> > >> D05_37_20957301_20141231.img.aux.xml
> >> > >> Size is 1000, 1000
> >> > >> Coordinate System is:
> >> > >> PROJCS["NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet_2011",
> >> > >> GEOGCS["GCS_NAD_1983_2011",
> >> > >> DATUM["NAD_1983_2011",
> >> > >> SPHEROID["GRS_1980",6378137.0,298.257222101]],
> >> > >> PRIMEM["Greenwich",0.0],
> >> > >> UNIT["Degree",0.0174532925199433]],
> >> > >> PROJECTION["Lambert_Conformal_Conic_2SP"],
> >> > >> PARAMETER["False_Easting",2000000.002616666],
> >> > >> PARAMETER["False_Northing",0.0],
> >> > >> PARAMETER["Central_Meridian",-79.0],
> >> > >> PARAMETER["Standard_Parallel_1",34.33333333333334],
> >> > >> PARAMETER["Standard_Parallel_2",36.16666666666666],
> >> > >> PARAMETER["Latitude_Of_Origin",33.75],
> >> > >> UNIT["Foot_US",0.30480060960121924],
> >> > >> VERTCS["NAVD_1988_Foot_US",
> >> > >> VDATUM["North_American_Vertical_Datum_1988"],
> >> > >> PARAMETER["Vertical_Shift",0.0],
> >> > >> PARAMETER["Direction",1.0],
> >> > >> UNIT["Foot_US",0.3048006096012192]]]
> >> > >>
> >> > >>
> >> > >>
> >> > >> In theory, these should all be EPSG:6543, so just assigning the
> >> correct
> >> > >> horizontal projection/datum with the epsg code should make things
> >> > >> usable. (i.e,
> >> > >> gdal_translate -a_"srs epsg:6543" --config GDAL_CACHEMAX 512 -of
> >> GTiff
> >> > >> -co COMPRESS=DEFLATE -co PREDICTOR=3 in.img out.tif ) ( Thank you
> >> for
> >> > >> the EPSG:6543 projection support in GDAL 2.0!)
> >> > >>
> >> > >> However, this only assigns the horizontal infomation, how does one
> >> > >> assign a vertical datum with a horizontal EPSG code?
> >> > >>
> >> > >> I see the NVD88 code in vertcs.csv , but how would I implement it
> in
> >> the
> >> > >> above command?
> >> > >>
> >> > >> Using gdal 2.0 Beta2
> >> > >>
> >> > >>
> >> > >> Doug
> >> > >> --
> >> > >> Doug Newcomb
> >> > >> USFWS
> >> > >> Raleigh, NC
> >> > >> 919-856-4520 ext. 14 doug_newcomb at fws.gov <mailto:
> >> doug_newcomb at fws.gov>
> >> > >>
> >> > >>
> >> >
> >>
> >>
> ---------------------------------------------------------------------------------------------------------
> >> > >> The opinions I express are my own and are not representative of the
> >> > >> official policy of the U.S.Fish and Wildlife Service or Dept. of
> the
> >> > >> Interior. Life is too short for undocumented, proprietary data
> >> formats.
> >> > >>
> >> > >>
> >> > >> _______________________________________________
> >> > >> gdal-dev mailing list
> >> > >> gdal-dev at lists.osgeo.org
> >> > >> http://lists.osgeo.org/mailman/listinfo/gdal-dev
> >> > >>
> >> > >>
> >> > >
> >> >
> >> >
> >> > --
> >> > Doug Newcomb
> >> > USFWS
> >> > Raleigh, NC
> >> > 919-856-4520 ext. 14 doug_newcomb at fws.gov
> >> >
> >>
> >>
> ---------------------------------------------------------------------------------------------------------
> >> > The opinions I express are my own and are not representative of the
> >> > official policy of the U.S.Fish and Wildlife Service or Dept. of the
> >> > Interior. Life is too short for undocumented, proprietary data
> >> formats.
> >> >
> >>
> >>
> >> --
> >> Spatialys - Geospatial professional services
> >> http://www.spatialys.com
> >>
> >
> >
> >
> > --
> > Doug Newcomb
> > USFWS
> > Raleigh, NC
> > 919-856-4520 ext. 14 doug_newcomb at fws.gov
> >
> >
> ---------------------------------------------------------------------------------------------------------
> > The opinions I express are my own and are not representative of the
> > official policy of the U.S.Fish and Wildlife Service or Dept. of the
> > Interior. Life is too short for undocumented, proprietary data formats.
> >
>
>
>
> --
> Doug Newcomb
> USFWS
> Raleigh, NC
> 919-856-4520 ext. 14 doug_newcomb at fws.gov
>
> ---------------------------------------------------------------------------------------------------------
> The opinions I express are my own and are not representative of the
> official policy of the U.S.Fish and Wildlife Service or Dept. of the
> Interior. Life is too short for undocumented, proprietary data formats.
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150609/323e8726/attachment.html
> >
>
> ------------------------------
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> End of gdal-dev Digest, Vol 133, Issue 17
> *****************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150609/9abaee1e/attachment-0001.html>
More information about the gdal-dev
mailing list