My goal is to check that a shapefile and raster have the same projection and datum using EPSG number.  One can easily extract EPSG number from a shapefile by the following:    ds=ogr.Open(r'....satates48.shp') utm_sr=ds.GetSpatialReference() utm_sr.GetAttrValue('AUTHORITY',1)   This works great but how does one do the same for a raster file, you cannot, to my knowledge get a spatial reference of a raster.  One can get projection, or projecitonref (I am not sure if there is any difference), but you cannot do more, like:  fn_raster=gdal.open(r'raster_file.tif') raster_sr=fn_raster.GetProjection().GetAttrValue('AUTHORITY',1)   This does not work. Perhaps there is an alternative way to do this?
120719
In Python, using GDAL, I've extracted a raster's projection as a WKT string as follows:  wkt = dataset.GetProjection() # wkt is 'PROJCS["GDA_1994_Transverse_Mercator",GEOGCS["GDA_1994",DATUM["GDA_1994",SPHEROID["GRS_1980",6378137,298.2572221010002],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]'   Using the WKT string, I can create a SpatialReference instance as follows:  src = osr.SpatialReference() src.ImportFromWkt(wkt)   This is easy-peasy. I can extract several parameters of the projection from src such as the UTM zone, etc., easily enough. But I can't figure out how to extract something like the name of the projection, i.e. "GDA_1994_Transverse_Mercator". This has surely got to be possible, but the Python API documentation may as well be nonexistent for all the use it is.  How do I extract the names of the projection and the geographic coordinate system?
100713
I have a HDF file that I want to convert to tif or ENVI file.  I know that the +Proj of the original file is: "+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m"  so:    gdal_translate -of "ENVI" -a_srs "+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m" HDF5 file.h5  output.img   but I found that the projection was not defined but instead I got this:         ENVI      samples = 1000      lines   = 500      bands   = 1      header offset = 0      file type = ENVI Standard      data type = 4      interleave = bsq      byte order = 0   My question how to deal with geo points?
104574
I try to create a common extent for rasters in QGIS 2.99. Here is my script:  import processing  processing.runalg("script:unifyextentandresolution",                   "/qgis_data/rasters/Image2.tif;                   /qgis_data/rasters/Image1.tif",                  -9999,"/qgis_data/rasters",True)   However, processing.runalg is being used for QGIS 2.99. I need it for QGIS 2.99
101291
When I use:  gdal_rasterize -b 1 -b 2 -b 3 -burn 255 -burn 255 -burn 255 polygon.shp tiff.tiff   I get a new tiff with a white hole in it with the shape of the polygon shp. However, when I use the reverse process:  gdal_rasterize -i -b 1 -b 2 -b 3 -burn 255 -burn 255 -burn 255 polygon.shp tiff.tiff   I get a blank tiff. I expect the old tiff shaped by the polygon (ie. raster values inside the polygon), instead.  Any ideas why this happens? is it possible to clip a raster with a shapefile in this way?
6090
I have Raster Data for a State which has Elevation Data in its pixels. I also have Vector layer which has the 300 divisions of State. what I need is, I want to get the average elevation value of particular division from the raster data using ArcGIS. I tried extract by mask in ArcGIS to extract individual polygon to clip the raster into 300 raster then from individual raster i have collected average elevation but it has taken days together to complete the work. Is there any other easy way to get it?
27250
I originally had an IMAGINE Land Cover raster file that I opened in ArcMap.   I used the Extract Feature from Raster to extract only the Developed class and saved it (still IMAGINE).   Then I used the Copy Raster to convert the IMAGINE to GeoTIFF, but when I try to view the GeoTiff in Google Earth or Photoshop, it doesn't display anything.   However, in ENVI, it displayed successfully.   What am I doing wrong?
80381
Could anyone provide with a way to calculate latitude and longitude from DEM, for each raster cell?   In previous posts about this, I have seen that some suggest calculating lat and long in certain sample point locations. This doesn't work for me, as I need complete rasters.   Also, I use ArcGIS 10.2 version so $$XMAP, $$YMAP does not work for me.   Is there any other way inside ArcGIS to get lat and long?
59487
I have a GeoTiff that has been supplied to me with non square pixels (I suspect due to reprojection to lat lon at some point). However the next step in my processing chain requires a Arc ASCII grid as input and so I need square pixels.   To save me trying to guess the original projection can any one tell me how to use GDAL to resample my image so the pixels become square?
102757
I use (1) ST_Length (ST_Transform (geom, 3395)) to get distances in meters from points in Geographic Reference Systems EPSG: 4326. The thing is I found that  with  (2) ST_Length (geom, true)  you can get also the distance in meters but using the spheroid. If I am right, is the option (2) more accurate if the GIS project use a Geographic Reference System for Global data, especially near to the poles? ST_Length reference:  EPSG:3395 reference:
122012
I have a huge amount of data in the form of CSV containing UTM coordinates as X and Y and an elevation value as Z information. I need to convert these data into a DEM as GeoTIFF for further analysis. In this case, a huge amount means 16 m. lines, with one point in X, Y and Z per line. The points are equally distributed, therefore an interpolation is not needed; each point just needs to be converted into a raster cell.  The original data came without separator, with fixed column widths. I already figured out how to convert the file syntax to use a separator instead of fixed widths and eliminate all space characters, using the stream text editor sed. From here on, normally my workflow would be to import the data into ArcGIS by creating a feature class from the X, Y and Z data and as a second step, converting the point shapefile into a GeoTIFF, using the Point to Raster tool. However, the file I currently have is way too big for this process.  Instead of the above described workflow, I was looking for an efficient alternative and discovered GDAL. However, in gdal_translate, the closest supported format I can find in the supported filetype list, is ASCII grid but no comma separated XYZ. Another difficulty is, that I have UTM coordinates, while most examples seem to use decimal degree coordinates. However, I need to stay within the UTM system (or at least, my output GeoTIFF needs to be in a UTM coordinate system).  So I am looking for a way to convert the CSV XYZ into a GeoTIFF, using GDAL, but so far wasn't able to find examples dealing with this exact problem. I would be very happy for some hints or even code examples.
55562
If I have, say, a rectangular and georeferenced raster file loaded in an ArcMap 10 document (.tiff w/ associated tfw), how do I easily find its centre point and store that point in a points vector layer?  Also, if I have multiple such rasters in my ArcMap document, how do I apply the process to all of them?  Unfortunately, I have zero Python experience. Therefore a programmatic solution is OK, but I will need specific instructions on how to load an existing script into ArcGIS 10, and run it on the rasters in questions. (BTW, the rasters are all in their separate layers).
105638
I have around 1,000 satellite images in tiff format, and I want to create a shapefile which will serve as an index to the rasters. This is something similar to a raster catalog, but I do not want to build a raster catalog.  Some obstacles I can forsee, is that the image are georefrenced, so they are not rectangular in shape (I am talking about the data area).  To clarify, I require the polygon to cover only the non-zero (or non nodata) pixels of the raster, and not the entire rectangular raster. Most of the answers so far, give a rectangular polygon, which covers the data, as well as the non-data pixels.  My Image     Result given by  tools I have examined (like raster catlog, various Arcscripts, custom Python script given in one of the answers):     Result that I want:
107836
I have 8 shapefiles in a folder with shapefile 1 and 8 in the same projection and the others (shapefile 2-7) if different projections.   Is there a way to use python to use shapefile 1 as the projection template to project shapefiles 2-7 in a new folder but also to only copy shapefiles 1 and 8 in the output folder without reprojection them but also not using hard values to keep the script flexible?  I have been trying to use batch project, but am having trouble excluding shapefiles 1 and 8, or only including shapefiles 2-7. Should I be using the Describe.SpatialReference to filter out?
76742
Shapefile: "table1.shp"   Field name: "Name" Field name: "Result"   Flat file: "referenceTable"   Field name: "STREET" Field name: "AKA"   I need to compare the value "Name" in "table1.shp" to the value "STREET" in "referenceTable"  If the values are the same then put the "AKA" value into the "Result" in "table1.shp"     Here is my code.  I got stuck on accessing the elements from the array  import arcpy, os, sys, string from array import *   table1 = "G:\\xStreetNew\\table1.shp" table1_column = ['Name', 'Result'] referenceTable = "G:\\xStreetNew\\Test.gdb\\referenceTable" referenceTable_column = ['STREET', 'AKA']  referenceTable_List = [] with arcpy.da.SearchCursor(referenceTable, referenceTable_column) as rCur: for r_row in rCur:     reference_row = r_row[0].split(",")     referenceTable_List.insert(0,reference_row) ##print sorted(referenceTable_List)  Table1_Update_Cursor = arcpy.da.UpdateCursor(table1, "Name") for row in Table1_Update_Cursor:   if row[0] == None:   row[0] = row[0] else:   row[0] = row[0].replace(";",",")   Table1_Update_Cursor.updateRow(row)  table1_List = [] with arcpy.da.SearchCursor(table1, table1_column) as tCur: for t_row in tCur:     row = t_row[0].split(",")     table1_List.insert(0,row) ##print sorted(table1_List)  Table1_row = len(table1_List) Table1_column = len([table1_List]) ref_row = len(referenceTable_List) for i in range(0,Table1_row):  #Table Row   for ii in range(0,Table1_column):   #Table Column     print table1_List[i][ii]   Below was my old post.  Sorry for the confusion.  I have a shapefile named "table1.shp".     I would like to get all values from fielname "Name" from this table and put them into a list. Able to output each element from the list     So far, I'm able to create a list but not too sure how to get to value "AVENUE 33" or "CHAPMAN ST"  This is my code:  import arcpy, os, sys, string from array import *   table1 = "G:\\xStreetNew\\table1.shp" referenceTable = "G:\\xStreetNew\\Test.gdb\\referenceTable" cursor = arcpy.da.UpdateCursor(table1, "Name") for row in cursor: if row[0] == None:     row[0] = row[0] else:     row[0] = row[0].replace(";",",") cursor.updateRow(row) table1_column = ['Name'] table1_curObj = arcpy.da.SearchCursor(table1, table1_column) table1_List = [] for row in table1_curObj:    table1_row = row[0].split(",")    table1_List.insert(0,table1_row)    print table1_List[0]
41217
I have a raster image (Sentinel-2 band 4) and a shapefile covering the same area than the image (same projection and extension). The shapefile has different polygons which are agricultural fields.   I want to create a new file where, for each polygon, the mean value of the pixels that are inside it is calculated.  Could you please help me to find an appropriate R code for that?
113293
I'm trying to write a script that can automatically use GDAL to crop the informational border from the FAA's published sectional charts (GeoTIFF) found :  I know that I can use gdal_translate -projwin ulx uly brx bry to crop it to a coordinate range, but I'm not sure where to get the values for the map area aside from by hand, or if there is possibly a better way to crop it.  Is the information necessary to crop the chart extractable from the files they provide?
34147
I'm trying to identify null/0 values in each of my raster bands. How can I loop through all the values of each band? I've tried this code, straight from the gdal/python cookbook website but this just gives stats for each band as a whole and for some reason it actually doesn't print a minimum value of 0 even though my raster does have 0 values, it shows a minimum value of 1.   I want to print every single cell value. Is there a way to do this?   from osgeo import gdal import sys  src_ds = gdal.Open( "INPUT.tif" )  print "[ RASTER BAND COUNT ]: ", src_ds.RasterCount for band in range( src_ds.RasterCount ):     band += 1     print "[ GETTING BAND ]: ", band     srcband = src_ds.GetRasterBand(band)      stats = srcband.GetStatistics( True, True )      print "[ STATS ] =  Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % ( \             stats[0], stats[1], stats[2], stats[3] )