1

I have a very large raster whose GDAL Info is as below:

Size is 172800, 67200
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,83.999167206000010)
Pixel Size = (0.002083333000001,-0.002083333000001)
Metadata:
AREA_OR_POINT=Area
ATTRIBUTE_LABEL=BDTICM_M
ATTRIBUTE_MEASUREMENT_RESOLUTION=1
ATTRIBUTE_TITLE=Absolute depth to bedrock (in cm)
ATTRIBUTE_UNITS_OF_MEASURE=cm
CITATION_ADDRESS=shanggv@bnu.edu.cn / tom.hengl@isric.org
CITATION_ORIGINATOR=College of Global Change and Earth System Science, Beijing Normal University / ISRIC - World Soil Information
CITATION_URL=http://dx.doi.org/10.1002/2016MS000686
CONFIDENCE_INTERVAL=M
DATA_FORMAT=Byte
DEPTH=None
DOWNLOAD_FTP_URL=ftp://ftp.soilgrids.org/data/recent/BDTICM_M_250m_ll.tif
HORIZON_LOWER_DEPTH=5000.00 m
HORIZON_UPPER_DEPTH=0.00 m
KEYWORD1=bedrock
KEYWORD2=depth
NO_DATA=255
PROJECT_URL=http://soilgrids.org
PUBLICATION_DATE=2017-03-10
RANGE_DOMAIN_MAXIMUM=5000
RANGE_DOMAIN_MINIMUM=0
SERIES_NAME=SoilGrids250m
SLD_NAME=soilgrids250m:BDTICM
SLD_URL=ftp://ftp.soilgrids.org/legends/BDTICM.sld
TECHNICAL_SPECIFICATIONS_URL=https://github.com/ISRICWorldSoil/SoilGrids250m
Image Structure Metadata:
COMPRESSION=DEFLATE
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-180.0000000, 83.9991672) (180d 0' 0.00"W, 83d59'57.00"N)
Lower Left (-180.0000000, -56.0008104) (180d 0' 0.00"W, 56d 0' 2.92"S)
Upper Right ( 179.9999424, 83.9991672) (179d59'59.79"E, 83d59'57.00"N)
Lower Right ( 179.9999424, -56.0008104) (179d59'59.79"E, 56d 0' 2.92"S)
Center ( -0.0000288, 13.9991784) ( 0d 0' 0.10"W, 13d59'57.04"N)
Band 1 Block=172800x1 Type=Int32, ColorInterp=Gray
NoData Value=-32768
Overviews: 86400x33600, 43200x16800, 21600x8400, 10800x4200, 5400x2100, 2700x1050, 1350x525`
I want to extract from it a part of data as xyz file. But I cannot seem to know how to do it. I have tried `gdal_translate -of XYZ -a_ullr -89.2972 60.0378 -117.719 44.901 "D:\\Mesh\\Shangguan et al. (2016) Bedrock data\\BDTICM_M_250m_ll.tif" "D:/Mesh/MRB/ascii shangguan for NCRB 0125 .xyz"

but it writes an xyz file more than 100GB and the size is increasing even more than 100GB as I write the question. Can I just have the data for "-89.2972 60.0378 to -117.719 44.901"?

  • You have asked about the same question with a bit other words gis.stackexchange.com/questions/245261/…. I consider this is a duplicate. – user30184 Jun 29 '17 at 7:16
  • @user30184 that question is if i can convert it to clip the same raster.. not for extracting xyz file..i need to do both actually – Kuljeet Keshav Jun 29 '17 at 7:20
  • Your original is compressed tiff tile COMPRESSION=DEFLATE and it holds 32 bit data Type=Int32. Such data as an uncompressed plain text format like XYZ takes evidently much more space. Your gdal_translate command does not make sense at all. Read gdal.org/gdal_translate.html -a_ullr ulx uly lrx lry: Assign/override the georeferenced bounds of the output file. What you need is -projwin which is documented on the same page. And think again what is upper left and what is lower right before giving values for -projwin. – user30184 Jun 29 '17 at 7:21
  • 1
    Did you already try gdal_translate -of XYZ -projwin -117.719 60.0378 -89.2972 44.901 "D:\\Mesh\\Shangguan et al. (2016) Bedrock data\\BDTICM_M_250m_ll.tif" "D:/Mesh/MRB/ascii shangguan for NCRB 0125 .xyz"? For comparison try also with GeoTIFF output. The output size should be so small that compression or tiling are not important. Remember from your other question that reading a block from a tiff file that is written row by row is inefficient and tiled tiff would be faster and take much less memory. – user30184 Jun 29 '17 at 8:17
  • @user30184 I tried gdal translate with projwin and it works! Though the file size is 4 gb but it gives the result. Thanks! – Kuljeet Keshav Jun 29 '17 at 18:26
3

XYZ is a verbose, uncompressed ascii format so you should anticipate large files. However, you have quoted your area as -a_ullr -89.2972 60.0378 -117.719 44.901. I'd have expected GDAL to give you an error message for those numbers since -117 is further west than -89 meaning you've actually quoted URLL and not ULLR. Perhaps because you are in WGS84, GDAL is content to go 'the long way round' the globe. I notice your prime meridian is Greenwich so perhaps you meant: -projwin -117.719 60.0378 -89.2972 44.901. This will give you an area of 26.52 x 15.13 degrees (instead of about 333.48 X 15.13! which would indeed be a massive file in xyz). Use -projwin to clip a subwindow. The -a_ullr parameter assigns new dimensions to your output but does not clip (so you may still expect a large file if you use that instead of projwin). For more information see the documentation here.

Even the area you intend in XYZ is going to be a big file especially with 15 decimal places for your X and Y values and int32 for your z values (though I'd not expect quite as much as 100GB+). Do you really need to use XYZ?

  • I need to find a way to get pixel values with their lat and long in a sheet so that I can process it with R or excel – Kuljeet Keshav Jun 29 '17 at 7:22
  • You can use XYZ, but just watch how you define the area as per the first paragraph in my reply. – MappaGnosis Jun 29 '17 at 11:07
  • Could you edit your answer a bit because -a_ullr just writes new metadata while -projwin is used for clipping? – user30184 Jun 29 '17 at 21:39
  • You are absolutely right, of course, and I should have read my reply back before pressing submit. Thank you for mentioning it. – MappaGnosis Jun 30 '17 at 6:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.