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