> > Message: 1 > Date: Tue, 9 Jun 2015 14:31:48 +0000 > From: "Ronquillo, Edgar Nahum" <[email protected]> > To: "[email protected]" <[email protected]> > Subject: [gdal-dev] Overlay shapefile onto Geotiff image > Message-ID: > < > 2086ce7e967c1b4f9d5897be1c3c746f31194...@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 <[email protected]>: > Send gdal-dev mailing list submissions to > [email protected] > > 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 > [email protected] > > You can reach the person managing the list at > [email protected] > > 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" <[email protected]> > To: "[email protected]" <[email protected]> > Subject: [gdal-dev] Overlay shapefile onto Geotiff image > Message-ID: > < > 2086ce7e967c1b4f9d5897be1c3c746f31194...@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" <[email protected]> > To: Even Rouault <[email protected]> > Cc: gdal dev <[email protected]>, Hermann Peifer > <[email protected]> > Subject: Re: [gdal-dev] assigning vertical datum to image files > Message-ID: > < > calqgvr2plpsskbeenswsk4qulus7t1wt9-g+thqhxq947-m...@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 <[email protected]> > 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 <[email protected] > > > > wrote: > > > >> Selon "Newcomb, Doug" <[email protected]>: > >> > >> > 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 <[email protected]> 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 [email protected] <mailto: > >> [email protected]> > >> > >> > >> > >> > >> > > >> > >> > --------------------------------------------------------------------------------------------------------- > >> > >> 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 > >> > >> [email protected] > >> > >> http://lists.osgeo.org/mailman/listinfo/gdal-dev > >> > >> > >> > >> > >> > > > >> > > >> > > >> > -- > >> > Doug Newcomb > >> > USFWS > >> > Raleigh, NC > >> > 919-856-4520 ext. 14 [email protected] > >> > > >> > >> > --------------------------------------------------------------------------------------------------------- > >> > 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 [email protected] > > > > > --------------------------------------------------------------------------------------------------------- > > 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 [email protected] > > --------------------------------------------------------------------------------------------------------- > 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 > [email protected] > http://lists.osgeo.org/mailman/listinfo/gdal-dev > > End of gdal-dev Digest, Vol 133, Issue 17 > ***************************************** >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
