texts: Operations between shape layer and csv layer

texts:
I have a .shp layer containing point geometry with xy coordinates in the attributes (left) and a .csv layer containing a square Origin/Destination matrix(right):

I need to perform this operation:

where:

x_i and y_i are the x coordinates of the origins

x_j and y_j are the x coordinates of the destinations

p_ij is the value in the O/D matrix.

Any idea on how to do this?


texts: What happens to .dbf file with TableToDBASE?

texts:
I have been trying to use TableToDBASE without success. I looked at this help link: http://desktop.arcgis.com/en/arcmap/10.3/tools/conversion-toolbox/table-to-dbase.htm. The code runs fine but it does not output anything. I can see in outLocation folder that something with the same name is being created, but when the code finishes running, the folder is empty. Here is my code:

import arcpy
workspace = r"C:\data\"
outLocation = r"C:\output"
arcpy.TableToDBASE_conversion("table", outLocation)

The output in IDLE is:

<Result 'C:\\Users\\Ariel Gomez\\Desktop\\todbf'>

What am I doing wrong? Why am I not getting a dbf file in the folder?


texts: Batch convertion of .tab files to a single .gml using ogr2ogr

texts:
My task is to convert .tab files from a directory (C:\Users\Vadim\Downloads\vec\MapInfo\snp) to a single GML file (snp.gml).
I use OSGeo4w Shell.
So I wrote a following command:

ogr2ogr -f "GML" C:\Users\Vadim\Downloads\vec\snp.gml -dsco FORMAT=TAB C:\Users\Vadim\Downloads\vec\MapInfo\snp

But it printed an error:

FAILURE: Unable to open datasource 'C:\Users\Vadim\Downloads\vec\MapInfo\snp' with the following drivers...

Could you please help me to get this command work.


texts: Converting MapInfo workspace/QGIS project to ArcGIS for Desktop?

texts:
I need to deliver maps in ArcGIS format to a client, but I only work in MapInfo and QGIS.

Is there any way to convert a MapInfo workspace or a QGIS project to ArcGIS for Desktop, keeping styles and thematic maps?

I don't have an ArcGIS for Desktop licence.


texts: Improving our spatial data by calculating "Close to Coastline", "Close to Park", "Close to Stadium" etc

texts:
A quick background for context is maybe helpful for this.
Essentially we have "vacation rentals" and I would like to add additional data to them automatically. We create booking websites for customers who have very niche sites, for example,

"Pet Friendly" rentals.
"Olympic Cottages"
"Surfing Adventure Lettings".

Quite often each of our customers always ask the same questions,

"It would be awesome if we could plot the distance and sort properties based on their relative distance from XYZ".

Usually XYZ is something like a list of dog friendly beaches, or cycle tracks, or wine vinyards. So esesntially points of interests which are easily definable. In those cases, we've always just asked the customer to provide a list of the lat/lngs of every possible "point of interest" that they would like us to calculate distance from".

More recently we've been asked to include "Close to Coastline".

What I want is to come up with a standard way of automating this. My ideas so far are to create a shapefile. For example, let's say we wanted to do "Close to Coastline", we could use something like

http://www.naturalearthdata.com/downloads/10m-physical-vectors/

At this point, I'm somewhat stuck in knowing how would be most efficient to calculate the distance from the edge of the coastline.

I'm thinking about being able to just create/add a shapefile per "concept" that the customer wants us to use for their product. But then do not know the process/tools available to me to know distance from each edge, or vector point.


texts:
I use ArcMap 10.5 and i have created a map with a layer with polygon-features. Each polygon has got its own name. I need all the labels to be shown in the Layout view. So far so good. My Problem is that some of the polygons are very small and some are big. It would be ideal if the labels adjust automatically to the size of the polygon and if that´s not possible i would like to "lock" the labels at a small size so that i can zoom in on them in an exported PDF-file. In the default-mode the label adjust to the zooming but if i zoom out some of the labels disappear.

How do I show all the labels in the finished map?


texts: A global, grid-like projection for creating heatmaps

texts:
I'm putting together an application where I need to create a vector grid which will be used to store and display a heatmap. It has the following requirements:

Can cover the entire planet.
The vast majority of grid squares won't have values.
I don't want to have to store the grid itself; I'd like to calculate it on the fly.
The scale of the data used with the grid could vary massively.
I anticipate wanting grid squares anything from 1km across to 100km across. (I'm aware how many this will be (~510million for 1km, ~51,000 for 100km)).
Values will be accumulated/aggregated for each grid square.
Ideally I'd be able to easily use smaller grid cells to calculate values for larger ones rather than store the larger grid cell values.
I'll be using OpenLayers to drape it over OpenStreetMap.
I'll be storing it in SpatiaLite or SQLite, so preferably supported by those natively (i.e. for SpatiaLite = a supported CRS; or for SQLite = a pure numbers based system).

So my question is: What projection should I use for this grid?

Also - is there a good way to designing this? Does anyone know of a good potential solution to this problem or have solved similar before? Or can point me in a useful direction.

Edit
Use-case - basically I'm aggregating up bounding boxes of various different shapes and sizes. They can be anything from a few hectares to thousands of square km in size. They may also be in different projections.

Below is a bespoke version of the sort of thing I'm aiming to accomplish automatically on a larger scale.

Many thanks.


texts: Transforming climate text data to raster in ArcGIS Desktop?

texts:
I have climate data in text files with no header in the format below from Udel's climate program.  Each row has a lat, lon, and a temperature for each month (14 columns total) from 1900-2010.  I need to generate a Raster layer for each year-month (1320 rasters) in ArcGIS Desktop 10.

I'm thinking the likely solution is a script to parse the existing text files to create an ESRI-formatted ASCII for each year-month and then use ArcToolbox to batch import ASCII to Raster.  However, I don't have any experience writing scripts.

Another option I considered is if there is a way to import a table or dbf to raster, but I couldn't find that operation.  Or, there may be another transformation that I'm not considering.


texts: Seeking QGIS equivalent of ArcGIS Erase tool?

texts:
Is there a QGIS equivalent of the Erase tool in ArcGIS Desktop?

The Erase tool:

  Creates a feature class by overlaying the Input Features with the
  polygons of the Erase Features. Only those portions of the input
  features falling outside the erase features outside boundaries are
  copied to the output feature class.


texts: Using ST_Intersection with three or more tables

texts:
I'm trying to achieve a multiple-table intersect with PostGIS while retaining the geom for each shape.

For any given property shape, I want individual cuts of that shape split apart by different characteristics of the below layers.

So for example:

480128 | 3Sh3W | Otamatea | MOZ | S | 2 | B' | 5 | P 5 | geom
480128 | 3Sh3W | Turangai| MOZ | L/S | 2 | B' | 5 | P 5 | geom
480128 | 3Sh3W | Makahu | ZOT | L/S | 2 | B' | 5 | P 5 | geom

The SQL I've used so far is:

CREATE TABLE sd_processed AS
SELECT
    A.title_no,
    B.erosion,
    C.series,
    C.domnzsc nzsc,
    C.ps,
    D.prd_class as prd,
    E.slope as slope,
    F.drain_clas as drain,
    G.veg,
    ST_DUMP(ST_INTERSECTION(
        A.geom,
        ST_INTERSECTION(A.geom,
            ST_INTERSECTION(B.geom,
                ST_INTERSECTION(C.geom,
                    ST_INTERSECTION(D.geom,
                        ST_INTERSECTION(E.geom,
                            ST_INTERSECTION(F.geom, G.geom)
                        )
                    )
                )
            )
        )
    )) geom
FROM
    sd_title A,
    sd_erosion B,
    sd_particle_size C,
    sd_potential_rooting_depth D,
    sd_slope E,
    sd_soil_drainage F,
    sd_vegetation G
WHERE
    ST_INTERSECTS(A.geom, B.geom) = true
    AND ST_INTERSECTS(A.geom, C.geom) = true
    AND ST_INTERSECTS(A.geom, D.geom) = true
    AND ST_INTERSECTS(A.geom, E.geom) = true
    AND ST_INTERSECTS(A.geom, F.geom) = true
    AND ST_INTERSECTS(A.geom, G.geom) = true
    AND ST_INTERSECTS(B.geom, C.geom) = true
    AND ST_INTERSECTS(B.geom, D.geom) = true
    AND ST_INTERSECTS(B.geom, E.geom) = true
    AND ST_INTERSECTS(B.geom, F.geom) = true
    AND ST_INTERSECTS(B.geom, G.geom) = true
    AND ST_INTERSECTS(C.geom, D.geom) = true
    AND ST_INTERSECTS(C.geom, E.geom) = true
    AND ST_INTERSECTS(C.geom, F.geom) = true
    AND ST_INTERSECTS(C.geom, G.geom) = true
    AND ST_INTERSECTS(D.geom, E.geom) = true
    AND ST_INTERSECTS(D.geom, F.geom) = true
    AND ST_INTERSECTS(D.geom, G.geom) = true
    AND ST_INTERSECTS(E.geom, F.geom) = true
    AND ST_INTERSECTS(E.geom, G.geom) = true
    AND ST_INTERSECTS(F.geom, G.geom) = true
LIMIT 1000;


texts: Taking classification image from ENVI and using it in ArcMap?

texts:
Loading ENVI file into ArcGIS Desktop? is a question similar to this but used a slightly different way of going about things so it made it not fit my situation so here's my situation.

I ran an isodata unsupervised classification. I have 8 relevant classes and all I need is to have all these classes editable in Arcmap. Unfortunately, the classification file is in an ENVI format and cannot be opened in Arcmap. I've looked at various answers to similar questions and they suggested using the ENVI Raster to Vector and then after that using a Vector to raster tool to make it editable in Arcmap. Which is fine and dandy but this image is huge and the process is incredibly time consuming. On iteration 4 of 8, a total of 2 hours of processing has passed and class 5 has twice as many points as 4.

So is there another way of doing this? Or is this going to be as fast as I can get this done? I've also toyed with the idea of using Arcmap's classification tool, but I don't know the processing on that.

If I am doing the most efficient method then I can wait... I would just like to know if there is a better way to do this before I have this process for 10 hours.


texts: Getting rid of DirtyAreas from topology using ArcPy?

texts:
I did a script the part of which creates and validate topology for a certain layer. Then topology object is deleted. But when it comes to edit session if I want to create new object it is shown that layer has dirty areas.
I tried to add lines like

arcpy.RepairGeometry_management(fc,"DELETE_NULL")
arcpy.AddSpatialIndex_management(fc)
arcpy.Integrate_management(fc,0.01)

But nothing works from script file or addin. The only solution is to launch tools listed above manually.
Is there any other way to remove dirty areas using code?

Here is the full code of my script.

import arcpy
import os

from arcpy import env
mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
arcpy.env.overwriteOutput=True
path = r"D:\Data\File.gdb\fd"
topologya = path + "\\Topol"
arcpy.CreateTopology_management(path,"Topol")
arcpy.AddFeatureClassToTopology_management(topologya,"lines")
arcpy.ValidateTopology_management(topologya)
arcpy.Delete_management(topologya)

arcpy.RefreshTOC()
arcpy.RefreshActiveView()


texts: OpenLayers: Modify features without virtual vertices?

texts:
When using the modify feature functionality of OpenLayers, I get the option to move the individual vertices around.
But I also get a set of virtual vertices - each vertex located on the edge between the original vertices. The virtual vertices allow me to kind of locally subdivide the existing polygon.

I can see that this feature is neat, but in many situations I don't want the user to be able to do local subdivision of the polygon - I want to maintain the original set of vertices while allowing the user to change the position of each of the vertices.
Is that possible with the existing OpenLayers controls?


texts: Convert polygon to raster and give value to null cells

texts:
Is there a way to convert polygon features to raster while also populating null cells with a specific value?  As in, can I create a raster from a polygon shapefile and have the tool populate all null cells with a value of zero (0)?  I'd like to do this because I will eventually sum several rasters using Raster Calculator, and it doesn't handle null values very well.

I know how to use the "Con(IsNull)" forumla in raster calculator to achieve this once the raster has been created, but would like to avoid having to perform this extra step if possible.


texts: Convert QGIS geometry into Shapely geometry to use orient method defined in Shapely

texts:
My idea is to fix the orientation of a polygon, so no matter the orientation it was created, when I run the code, the polygon become clockwise. How is it possible?

I searched and found the shapely library, and this documentation, may help, but I don't know how to implement.

layer = iface.activeLayer()
features = layer.getFeatures()

for feature in features:
    shapely.geometry.polygon.orient(polygon, sign=1.0) # How to implement?


texts:
I have a small problem with the Python console in Qgis. I am working through some code .py but lately the console does not show results, for example I write directly on the cosola

print "prueba"

Enter and it shows nothing, what could be the problem?


texts: Converting Geo-normalized SVG data to geojson

texts:
I'm working on a project for which our team geonormalized some TIFF files with MapWarper, then used photoshop to vectorize only the buildings within the images. We now have svg data for which we know the geo-coordinates of the bounding box (that is, the geonormalized representation of the original image corners).

My question is: Is there a way to convert an svg file for which one knows the coordinates of the bounding box to geo json? Any help others can provide on this question will be very welcome!


texts: How to configure "-nlt MULTILINESTRING" for DXF export?

texts:
When I try and export a polygon shapefile to DXF, the resulting features in AutoCAD are hatch patterns, which is as expected if you read http://www.gdal.org/drv_dxf.html, "Polygon and MultiPolygon features are written as HATCH entities".

After doing some more digging I found that OGR2OGR can export a polygon shapefile to polylines using the following:

"C:\OSGeo4W64\bin\ogr2ogr.exe" -f "DXF" some.dxf some_shapefile.shp -nlt MULTILINESTRING

Can I somehow add "-nlt MULTILINESTRING" to the Custom Options boxes using the "Save As" function in Qgis?


texts: How to update annotations from fields values to annotations objects

texts:
I have annotations feature class and need to update some field offset and alignment values

When I edits them manually at attributes table in ArcMap, new values applies to objects and they changes thrir positions at map.

And when I updates them via script:

with arcpy.da.UpdateCursor(annotations_fc, ['VerticalAlignment', 'YOffset']) as cursor:
    for row in cursor:
        row[0] = 3
        row[1] = 2.0
        cursor.updateRow(row)

del cursor

attributes table has new values but these values doesn't apply to map objects

Arcpy has UpdateAnnotation_management method for annotations updating. But as I understand, it updates field values from annotation objects when I needs reverse action - update annotation objects from fields values.

Is such method exist? If not in arcpy, at least in arcobjects.net?

Or am I doing something wrong?


texts: Why does this python function return a value of 0 in field calculator?

texts:
I am trying to calculate the percentage of a total which the median value makes. I have a column with the median column labelled (i.e. 1st,2nd,3rd,4th column is median), the columns of values and the total of all values. I wrote a small python function in the Codeblock and it returns a value of 0 rather than a percentage. All the fields involved are dobule. I know I did this correctly before but I've tried almost everything to fix my code and cannot get it to return anything other than 0. Here is a screenshot of what I did in field calculator:

And here is the code written out:

def medfind(peri1,peri2,peri3,peri4,peri5,peri6,peri7,index,tot):
      dict = {1:peri1,2:peri2,3:peri3,4:peri4,5:peri5,6:peri6,7:peri7}
      return dict[index]/tot

I have double checked all the fields and everything does have a value. I just can't figure out why it's not working. This is in a .shp by the way, and I have previously used this function to success.


texts:
I have created a buffer layer and I have a points layer.

How do I count how many points lie in each buffer using ArcMap?

I have been using a spatial join and I know I should probably use the Sum field feature but I don't really understand how that works and if that's the correct solution.


texts: Using Google Earth imagery data for portable solution?

texts:
I wish to clarify with users who are using Google Earth portable solution. As per portable solution link http://www.google.com/enterprise/earthmaps/portable.html, we can create part of globe and create offline map to use it.

However FAQ section in : http://www.google.com/enterprise/earthmaps/earth_enterprise_faq.html mentioned as below:
Does the Google Earth Portable solution allow you to use Google's imagery?
No, the Google Earth Enterprise portable solution is intended to be used in situations where your employees need to access your internal imagery and mapping data. Google does not provide customers with Google licensed data for portable use. Google has many partners that provide custom Google Earth Enterprise portable solutions with high resolution imagery.

Then i am not sure how to use portable solution in offline without Google's imagery data. Any experiences or details on this solution are welcome.


texts: Read N number of rows from shapefile using GeoPandas

texts:
Pandas has a parameter nrows for selecting how many rows will be read into the DataFrame.
Does GeoPandas have a similar parameter? Couldn't find in the documentation.


texts:
How to convert kml to shapefile using python code? I'm using the code below but it's not properly coverting:

vlayer = QgsVectorLayer("D:\KML-Images\dp0844061p.kml", "line", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)
QgsMapLayerRegistry.instance().removeMapLayer(vlayer)

_writer = QgsVectorFileWriter.writeAsVectorFormat\
(_vlayer,"hoppla.shp","utf-8",None,"ESRI Shapefile")

_writer = QgsVectorFileWriter.writeAsVectorFormat(vlayer,r"D:\KML-Images\dp0844061p.shp","utf-8",None,"ESRI Shapefile")

_writer = QgsVectorFileWriter.writeAsVectorFormat(vlayer,r"D:\KML-Images\ShapeFile\dp0844061p.shp","utf-8",None,"ESRI Shapefile")


texts:
I am using QGIS 2.18.3 to create an atlas. This works well and I now want to call some functions to create richer content in a text box. To do this I am calling a python function from an expression. Works fine.

Since I am working on several computers I would like to keep track of my python code. I am using GitHub for this. However - QGIS takes a global approach and by default places my code in this folder (yes using Windows):

C:\Users\raglar\.qgis2\python\expressions

I would like to keep the python code related to this project relative to this particular QGIS project file. This way I will be able to keep track of the functions relevant to this particular project. Is this possible?


texts: Measuring distances between many points in ArcGIS for Desktop?

texts:
So I'm new to this, the answer is probably really simple.

I have a point map of X Y data points (one for each participant) for a certain activity.  I also have a shape file of zip codes, which I then made a point file out of this with the centroid for each zip code.  I want to measure the distance from the point for the activity to the zip code that corresponds to each specific participant.  I joined the attribute tables for these so that each participant has the a column for the lat, and long for the activity and the zip code lat and long.  I must be close.

How do I ask ArcGIS to measure the distance between these points?

I'm using ArcMap 10.2.


texts: Automated creation of perpendicular lines between a point layer and a line layer

texts:
I am using QGIS and I am looking for a script, or a plugin, that can create a large number of perpendicular lines from a point layer to lines in a separate line layer.

So far, I've tried to use the Hub Distance function in MMQGIS (turning the lines into points and then connecting the points to the nearest hub) and the respective tool of the QGIS geoalgorithms. Neither worked. Both take more than 2 hours and create either lines all of the layer or lines that are not perpendicular or connected to the points.

In the picture, you can see the present status of the project. The perpendicular lines should run from the points to the nearest line. In the end, I would like to use intersection points with a line between the points and the country borders to create a buffer of 4-sided polygons that is two polygons deep. I mention this in case that there is an easier way of doing this.

I know that there are some posts on how to create perpendicular lines, but neither of them solved my problem.


texts: Trying to extract a list of Unique Values from a field using python

texts:
I have a number of columns in a number of tables withinh a FGDB where I need to extract the unique values for each column.

For Example: the values may be  [1,2,2,2,3,4], and I am trying to return [1,2,3,4]

I could do this job a number of other ways in ARCGIS  but I am trying to extend myself.

I have found a piece of python on the web that I think will do the job but I am struggling to get it to run (I keep getting an invalid syntax error as I keep getting the syntax error in line 3)  this will no doubt be a really simple user error.

Code Snippet below

import arcpy

def unique_values(r'N:\GISProjects\Landuse\Plant_Biosecurity_Project\ArcGIS_Online.gdb\Holdings_Property_Merge' , 'LU_ALUMMaj'):
with arcpy.da.SearchCursor(table, [field]) as cursor:
    return sorted({row[0] for row in cursor})

This is the Error Message I get from sublime text:

 File "C:\Users\hawkinle\Desktop\STDTAS\Unique_Data.py", line 3
def unique_values(r'N:\GISProjects\Steve_Eastwood_Landuse\Plant_Biosecurity_Project\ArcGIS_Online.gdb\Holdings_Property_Merge' , 'LU_ALUMMaj'):
                                                                                ^
SyntaxError: invalid syntax
[Finished in 0.1s with exit code 1]

Updates from Original Question

I have now updated my code with the answer provided by below but am recieving a secondary error.

New Code Snippet:

import arcpy

def unique_values(table , field):
    with arcpy.da.SearchCursor(table, [field]) as cursor:
        return sorted({row[0] for row in cursor})

myValues = unique_values(r'N:\\GISProjects\\Landuse\\Plant_Biosecurity_Project\\ArcGIS_Online.gdb\\Holdings_Property_Merge' , 'LU_ALUMMaj')

print (myValues)

I am getting a new error message related to a runtime error

Traceback (most recent call last):
File "C:\Users\hawkinle\Desktop\STDTAS\Unique_Data.py", line 7, in <module>
myValues = unique_values(r'N:\\GISProjects\\Steve_Eastwood_Landuse \Plant_Biosecurity_Project\\ArcGIS_Online.gdb\\Holdings_Property_Merge' , 'LU_ALUMMaj')
 File "C:\Users\hawkinle\Desktop\STDTAS\Unique_Data.py", line 4, in unique_values
with arcpy.da.SearchCursor(table, [field]) as cursor:
RuntimeError: cannot open 'N:\\GISProjects\\Steve_Eastwood_Landuse\\Plant_Biosecurity_Project\\ArcGIS_Online.gdb\\Holdings_Property_Merge'

[Finished in 8.0s with exit code 1]

I assume from the reading I have done that this relates to setting the env.workspace ?


texts: OGR - 'NoneType' object has no attribute 'CreateLayer'

texts:
I have an error when i wan't to create a Shapefile layer with OGR. When I run my script for the first time, everything all right, layer is well created. But when i run it for the second time CreateDataSource() is None and CreateLayer() returns an error. If i run it a third time, it works again.
If i look at the documentation everything seems ok. How to explain and prevent this ?
import osgeo.ogr as ogr
import os

output="full/path/to/my/layer.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(output):
    driver.DeleteDataSource(output)
ds = driver.CreateDataSource(output)
layer = ds.CreateLayer(output, geom_type=ogr.wkbPolygon)


texts: Writing shapefile metadata with geotools

texts:
Using GeoTools, I'm able to write features to a shapefile by essentially passing a SimpleFeatureCollection to a ShapefileDataStore. The feature collection provides both a schema and features including their attribute values.

In addition to having the row-per-feature attribute table in the shapefile, I'd like to include file-level metadata - either structured fields or at least free-text title/description fields. Is there support for doing this in GeoTools? Or do I need to serialise my own *.shp.xml file?


texts: Understanding group attribute in mapserver's cluster?

texts:
I am trying to use the Group function in Mapserver's Cluster.

We have a cluster defined where when we have a lot of objects on top of each other we group them in a cluster.  Show a circle and a count.  However the cluster is made up of two types of locations.  And we would like to see two numbers (a number for each type) in the cluster representation.  We were trying to use Group for this but have not had any luck and cannot find an example.


texts:
I have obtained LiDAR data from my environmental agency in order to obtain the relative heights of each watercourse station. A far as I know, LiDAR signals do not get returned over water bodies. However my data does in fact show many last returns. Would these returns be reflections from nearby objects (vegetation)? The height values do however show that there is a relative drop from the river bank returns to the points returns on the water surface.

The image shows the course of a river and the LiDAR surface returns.

Could anyone explain why I get returns over water bodies and whether these returns have any practical use?


texts: cropping gtif with gdal_translate or qgis clipper

texts:
I have a Gtiff which loads into QGIS and seems to be OK.  It's a US topo.  When I try to clip it to remove the margins it fails (from QGIS or from cmd line) with this error:

dt1:gdal_translate -projwin 253.875 38.625 254 38.5 -of GTiff  CO_Salida_West_451198_1983_24000_geo.tif topo.tif

Input file size is 6579, 8096
Computed -srcwin 3441 7202 0 0 from projected window.
Error: Computed -srcwin 3441 7202 0 0 has negative width and/or height.

I also tried it with the corresponding negative value for longitude
The output from gdalinfo:

Driver: GTiff/GeoTIFF
Files: CO_Salida_West_451198_1983_24000_geo.tif
Size is 6579, 8096
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.9786982139006,
                AUTHORITY["EPSG","7008"]],
            AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4267"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",33],
    PARAMETER["standard_parallel_2",45],
    PARAMETER["latitude_of_origin",38.5],
    PARAMETER["central_meridian",-106.0629999999992],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-6739.200995264532139,14674.596086209008718)
Pixel Size = (2.032000000000000,-2.032000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (   -6739.201,   14674.596) (106d 8'26.95"W, 38d37'58.41"N)
Lower Left  (   -6739.201,   -1776.476) (106d 8'26.36"W, 38d29' 1.98"N)
Upper Right (    6629.327,   14674.596) (105d59'11.22"W, 38d37'58.41"N)
Lower Right (    6629.327,   -1776.476) (105d59'11.80"W, 38d29' 1.98"N)
Center      (     -54.937,    6449.060) (106d 3'49.08"W, 38d33'30.29"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 3290x4048, 1645x2024, 823x1012, 412x506, 206x253, 103x127, 52x64
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 3290x4048, 1645x2024, 823x1012, 412x506, 206x253, 103x127, 52x64
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 3290x4048, 1645x2024, 823x1012, 412x506, 206x253, 103x127, 52x64


texts: How to use non-ESRI tiles inside Arc Online

texts:
This link now appears to be outdated.  LINK

Is there another way to display other tiles as basemaps with the current configuration?

Here's the link I'm trying to use:  http://b.tile.stamen.com/watercolor/12/655/1582.jpg


texts: estimate travelling time for the route using pgRouting and OSM data

texts:
Could someone share with me some ideas or your experience how to estimate travelling time based on the provided route within pgRouting? The problem is that roads inside and outside cities have the same type (for example "highway=primary") but different max speed allowed. But this max speed attribute is not always set in real OSM data.
The only possible solution I see at the moment is to preprocess all roads before calculations and set for each some expected speed. For example, if "highway=primary" and it's inside the city then expected speed is 40km/h, if it's outside then expected speed = 90 km/h. But how to detect that the road is inside some city bounds or not?
Any ideas are welcomed!


texts: Convert UK OS grid coordinates to decimal degree WGS84?

texts: I've got a list of coordinates in the UK Ordnance Survey grid system e.g. "NZ 27 41" or "NZ 284 449". How do I convert them in QGIS 2.10 to decimal degrees so I can map them in QGIS in WGS84? BTW, is there a way to deal with the spaces in those coordinates?


texts: Find optimal chunk size using opt_chunk_size() in lidR catalog?

texts:
Is there a way to find the optimal chunk size using opt_chunk_size() when using a lidR catalog?

I seem to always randomly choose integer values until the program quits throwing errors at me. For example, this is from the catalog_apply() lidR documentation (p.11):

## ===================================================
## Example 2: compute a rumple index on surface points
## ===================================================
rumple_index_surface = function(cluster, res)
{
las = readLAS(cluster)
if (is.empty(las)) return(NULL)
las <- lasfiltersurfacepoints(las, 1)
rumple <- grid_metrics(las, ~rumple_index(X,Y,Z), res)
bbox <- raster::extent(cluster)
rumple <- raster::crop(rumple, bbox)
return(rumple)
}
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
project <- readLAScatalog(LASfile)
opt_chunk_buffer(project) <- 1
opt_chunk_size(project) <- 120 # small because this is a dummy example.
opt_select(project) <- "xyz" # read only the coordinates.
opt <- list(raster_alignment = 20) # catalog_apply will adjust the chunks if required
output <- catalog_apply(project, rumple_index_surface, res = 20, .options = opt)
output <- do.call(raster::merge, output)
plot(output, col = height.colors(50))


texts: How to create a raster (topo to raster) with non-rectangular shape?

texts: I'm using Arcmap 10, I need to create a raster from contour lines. These contour lines were clipped out with a certain shape with the purpose to modify them only within that shape and then create the raster from these modified contour lines. The problem is that when I come to the TopoToRaster step, the resulting raster is rectangular, while I'd need it with the shape I used for the clipping. How can I do this? Should I change some options in the Spatial Analyst? How?


texts: Importing recent NAVTEQ/HERE RDF Data into PostGIS?

texts:
I am trying to import a recent version of navteq data into a postgres DBMS. I have downloaded the 2014Q2 version of the data from https://here.flexnetoperations.com as well as Europe_RDF_rdf_customer_software. The RDF Software has a java program to import it into Oracle. I changed the JDBC settings to point to the Postgres DBMS. This did not help a lot, since the script uses sqlplus, a component of Oracle.

Are there any tools (from navteq/here) that could import all the navteq RDF data into postgres?


texts: Using Linkage mapper to build corridor between core forests?

texts:
I build corridors between core forests using the Linkage mapper (https://circuitscape.org/linkagemapper/). For inputs i used core forest (Polygon file) and resistance raster as 1 (uniform resistance with value 1). Euclidean tool was used for network adjacency method. I also limited the maximum distance between core forests for corridors as 20 Km meaning forests farther than 20 Km will have no corrdiors. With these specifications, the tool was run and the results for corridors were obtained.

But i have difficulty understanding the final raster values of the corridors obtained in corridor.gdb.

What are those values in the raster file and how did the tool reach to the result?

I am attaching the screen shot of the output file with corridor raster along with LCP(least cost path marked as red).

The values of corridors are bit strange and why is it so?

Could someone help me in explaining the image?

The LM2903_corrdiors is the final output image with values called  the normalised least cost corrdior value which is equal to Cost weighted distance from core forest (a) + Cost weighted distance from core forest (b)-Least cost distance when moved along the ideal path. (I made the symbology such that values greater than 5000 m gets no colour). I found this equation from the guidelines which i didn't understand completely. The marked red line is the least cost path and takes a strange shape even when the resistance is uniform.

Why is it so?


texts: Get Output Preview from QGIS ExpressionBuilder Dialog

texts:
I am trying to determine the data type of a field dynamically using the output preview (highlighted in the image below). However, I couldn't get a method in QGIS 2.18 API that does that.

What I am looking for is the output preview value when the preview is shown.


texts: Accessing feature attachments using ArcObjects?

texts:
I am attempting to access a feature's attachments in ArcObjects and I'm following the docs for Working with feature attachments.  Everything is fine until attempt to access IAttachment.Data because I am expecting there to be data, but that property is null:

which means I can't proceed past:
IMemoryBlobStream memoryBlobStream = attachment.Data;

I am thinking there is a missing step (or steps) in the docs that would fill the attachment object's Data property with actual data? How is this done?

I am able to access the attachment via the Identify window within ArcMap, so I know the data is there:


texts: create convex hull polygons / cluster for sets of N points in a point dataset

texts: I have 50k points in a feature class which I need to build convex hull or polygon for each set of 50 points that are nearby each other. I guess the question is very similar to create census areas for each 200 households. I've tried all tools in ArcGIS Desktop 10.1 but none of them does what I'm after. I also tried some other third-party tools which allows me to create convex hull based on K-means algorithm - means I can specify the number of convex hulls to be generated but still can't see how to create convex hulls for a specified number of points. I'm open to any solutions not just limited to ArcGIS, but the final result need to be imported back to ArcGIS. Any advise would be appreciated! Thanks.


texts: PyQGIS: How to get field values from a vector layer based on an array of field names

texts:
I have a list with the field names and I am getting the values from the vector layer

shp = "D:\file.shp"
from osgeo import osr,ogr, gdal
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(shp,0)
layer=dataSource.GetLayer()
list=['field1','field2','field3']

I am getting the follow error if I am doing it in the list

ValList = [feature.GetField(j) for j in list]

raise ValueError("Illegal field requested in GetField()")
ValueError: Illegal field requested in GetField()

So what I have to do if I want to get the values in the list?


texts: Convert CAD to graphics for layout template in ArcGIS

texts:
I've received a CAD template and need to use this as an ArcGIS layout template. I can go to layout view in ArcGIS and "convert features to graphics" as suggested in this answer.

This answer suggests redrawing all the lines in the CAD template (step 5). Converting CAD features to graphics still means that the graphics have a "spatial" location and will be visible in the data view. Is there a way to divorce a CAD conversion from any spatial information and effectively have it behave as lines drawn using the "draw" toolbar in the layout view?


texts: Getting different X coordinates for points which are displayed right next to each other?

texts:
I want to connect two points with a line with the "XY to line" tool. Therefore I add the XY coordinates with the "add XY coordinates" tool. Although the points are right next to each other the X values of one of the point layers are far too high. I can tell which layer is causing the problem as when I run the "XY to line" tool the line starts at the right point but goes to the middle of nowhere.

See the images below.

I already fixed it, I projected one of the point layers to the wrong coordinate system.


texts: Snapping to vertices in vector tile layer using ArcGIS Pro?

texts:
I am trying to create features in a geodatabase layer using the tracing tool of ArcGIS Pro. I am able to use other shapefile layers in the project for tracing but vector tile layers loaded from the portal are ignored by the tool.

Is it possible to snap to a vector tile layer?


texts: Error when running GeoServer with GWC S3 BlobStore plugin installed

texts:
I have GeoServer 2.13 installed and I would like to use the GWC S3 BlobStore plugin which I downloaded from here:

http://ares.opengeo.org/geoserver/master/community-latest/geoserver-2.13-SNAPSHOT-gwc-s3-plugin.zip

I followed the instructions from here:

http://docs.geoserver.org/latest/en/user/community/gwc-s3/index.html

I unzipped the contents into WEB-INF/lib directory as per the instructions however when I startup GeoServer I get a bunch of errors, here is the log file:

https://pastebin.com/BYdR4jp9

I believe this question has been asked a few years ago, but I couldn't find it. I remember the answer being that the nightly build was problematic. Does anyone know if this is the case again or how to fix the problem?


texts: Integrating Google Maps Traffic and Transit in QGIS

texts: Is there some manner to download and shown Google Traffic and Transit layers in QGIS? Something made with a Plugin or inserting a shapefile.


texts: Static traffic visualization for certain streets with the help of historical traffic data, in QGIS?

texts:
The goal is to make a static map highlighting certain streets in a city that usually have high traffic (for example, highways and main roads).

I am using OpenStreetMap for the city. The streets need to be highlighted by various colours, depending on their average yearly traffic.
For instance,

green colour for low traffic streets
orange for medium traffic, and
red for high traffic streets.

I have the data of average yearly traffic for these streets, which I have classified into low, medium and high traffic segments, in MS Excel.

I am new to visualization. I tried to search this online on various platforms but have no luck as of yet. Any suggestions to what I might try/look at will be useful. I am using QGIS 3.2.


texts: ArcObjects: Is there a way to insert an error into the Start Editing dialog

texts:
When I use the the Editor Start Editing toolbar item, a dialog with information, warning and error appears before the edit session is started:

Is it possible to insert programatically my own error in the dialog, so that it is shown and prevents he user to start the edit session? (Description)


texts: Free ServerContext has failed because the Wait Timeout Interval has elapsed

texts:
ArcGIS Server 10 SP4 logs below entry in log file. We  have a single map service with 100+ layers which has Maximum time client can use a service is 600 seconds (10 min).
I am really surpised when some of requests failed with below message. I could not understand (or believe) which request to map server is taking more than 10 min. Can some one throw some light on this? I read all ESRI documentation and forums.

A request to obtain a free ServerContext has failed because the Wait Timeout Interval has elapsed.


texts: Joining multiple polygons inside a single polygon in QGIS

texts:
I need to somehow "join" information from two shape files:

zoning - these are large polygons indicating the zoning code for a moderately large area, typically hundreds of separately owned parcels (see #2) are fully contained in each. The boundary files of zoning seem pretty "rough hand drawn", with boundaries typically running down the MIDDLE of the streets. They are however completely closed and do not overlap.
parcels - these are relatively small polygons, with precise lot boundaries, one for each separately owned parcel. None of these polygons extend out into the street. Related to each parcel id is a lot of useful tabular data, however.... that data doesn't include the actual zoning code (eg residential, commercial, industrial, office/residential, etc, etc, etc.) Ideally I'd like to just "compute" or "look that up" from the more aggregate zoning shapefiles and add a column to the parcels tabular data (unless of course there is a better way...)

Totally new to vector->geo processing tools->

Gave Union a quick whirl but got odd results, long processing time, and a very large output file, perhaps due to the fact the darn zone file borders run down the center of the street vs. the strict property boundaries contained in the parcels file.


texts: How to use Orthoengine in PCI Geomatica to georeference Pleiades 1B TIFF to WGS 84

texts:
I am trying to use the orthoengine in PCI Geomatica 2015 to reference an image to WGS 84 (ESPG: 4326).

Acquisition mode PX
Spectral processing PA+XS (resolution is 70cm)

I have the associated RPC file.

The problem I have is that when I try to set the projection to that in the RPC file (WGS 84) I get an error message stating that the output georeferencing must be projectable.

Any suggestions?


texts: ODBC connection issue on 64bit architecture

texts: I'm migrating an addin to a 64bit machine using Arcmap 9.3.1. Upon execution, it returns the following error - "[Microsoft][ODBC Driver Manager] The specified DSN contains an architecture mismatch between the Driver and Application". I'm aware through sleuthing that there exists two versions of the odbcad32.exe; one in the System32 folder, and one in sysWOW64. Oracle 11g connections are already set up such that one should be able to configure either driver in the control panel, but upon attempting to configure the necessary driver using the 32 bit version, the very same error occurs. Only the 64 bit tests correctly, but the software uses 32. My question is then - is it possible to point to odbcad32 in the system32 folder, and otherwise, why is it that the ODBC driver might be returning the architecture mismatch error upon configuration in the control panel? I also tried rebuilding the tool source (VB .NET if that matters) in x64 but this did nothing.


texts: Setting the initial visibility of the KML NetworkLink accessing ArcGIS Server map service

texts:
The objective is to create a KMZ for Google Earth that I can serve up that accesses an ArcGIS Server map service with KML capability turned on (KMLServer) and control the initial visibility of the layers. However, NetworkLink does not honor the ArcMap visibility setting when set in the original map document (mxd); everything is simple set to visible on initial load.

Anyone have any idea(s) how I can control the visibility of the NetworkLink's sublayers (on initial load)?

I've considered using multiple NetworkLinks (as seen below) and using <visibility> but that solution isn't ideal since it creates a bit of a nested mess.


texts: Getting result in metric from QGIS Distance Matrix tool in QGIS?

texts:
I'm trying to calculate the distance between points in layer with distance matrix tool. And keep getting strange results no matter what. Like this
MASTER  ------------SLAVE  ---------------- DISTANCE  --------- SHOULD BE
Metrotown ---------Knox Mountain ------ 0.036069979------------ 48.2 km
Metrotown -------- Central  --------------- 0.048071058------------- 29.3 km
Metrotown -------- Rosedale ------------ 0.050883421-------------   60.1 km
Metrotown -------- Kelowna -------------- 0.051697257------------   340.4 km
Metrotown -------- Ocean Park --------- 0.054471001-------------    29.1 km
Metrotown --------  Langley --------------  0.055017692-------------    18.7km

 First my thought was that the distance is in degrees not in meters. So I tried different CMSs but it didn't change anything. Than I tried to duplicate layer and calculate distances between points between two layers. But nothing changed.
NAD83 CMS gave me another unrealistic result ( smth like 5665565, 5675656, etc)

What CMS I should use?

I even tried to do this in postgis with ST_Distance and still getting same results.


texts: How to group by attributes and give unique categories in QGIS?

texts:
I have features and each of them has an attribute about what it is. In my case, I have points representing trees and each tree has its species name as an attribute.

I would like to group them by species name and using numbers from 1 to n as seen on picture and if possible, not by hand because sometimes I have hundreds of species. I was thinking of some kind of expression which includes $rownum and tried a few things but no results. SQL queries also came to mind, but I have limited experience with them and I've never given them a try in QGIS.

I really don't want to export data to a spreadsheet, I would like to solve this in QGIS if possible.


texts:
My question is an extension of Vertical lines in a polygon shapefile. Kindly refer to that question first.

What you will see is a method of generating vertical lines with respect to the bounding box, at user-defined spacing. I understand that OGR, Fiona, Shapely etc. can be used to do the next step of clipping, but I do not understand their utilization.

How do I read one line of a polygon shapefile? Every application that uses Shapely shows how to generate the LineString, Point or Polygon but never to read an existing shapefile

Kindly assist me with at-least a skeleton structure so I can build on it.


texts: Making map where overlapping symbols increase Opacity?

texts:
I was looking at some maps like: Map from Ship Logs

or the Facebook worldwide friendships Mapped (via the Flowing Data blog)

These Maps look good because single lines are almost transparent, and as more and more lines are drawn, the opacity goes up, and you get patterns only with lines where multiple lines are present.

I have some data of Bus Routes in shapefile format, and would like to do something similar from this.

Which software will allow me to do this? Which software can create this effect?


texts: Generating a neighborhood pixel raster from existing raster using PostGIS

texts:
I am starting with a PostGIS raster table containing anywhere from 6-12 bands.  In a nutshell, the goal is to take each of these bands, perform a 3x3 neighborhood average on each pixel in the band, and output a new 12-band raster with the neighborhood-averaged bands and the original geometry intact.  I could likely do this with pure GDAL, but the rasters at this stage in the process have already been imported into the database and the original GeoTIFF files deleted, so using GDAL would likely require an infrastructure change.

My current concept uses ST_MapAlgebra to take the bands from the source raster in one-by-one and perform the 3x3 average to return the single averaged band.  I then would use ST_Union on a setof these averaged bands to create a single raster.  I then need to somehow generate a new table from the output of the ST_Union.  I am trying to do this in as few steps as possible, so ideally the process would take 1-2 queries.

My ST_MapAlgebra callback function is this:

CREATE OR REPLACE FUNCTION
                    callback_fn(pixel float[], pos integer[], variadic userargs text[])
                        returns float
                        language plpgsql
                        immutable
                    as $$
                        declare
                            avg float;
                            x integer;
                            y integer;
                            p integer[];
                        begin
                            x = pos[1][1];
                            y = pos[1][2];
                            p = pixel[1];
                            avg := (p[x-1][y-1] + p[x][y-1] + p[x+1][y-1]
                                    + p[x-1][y] + p[x][y] + p[x+1][y]
                                    + p[x-1][y+1] + p[x][y+1] + p[x+1][y+1]) / 9.0;

                            return avg;
                        end
                    $$;

My current attempt at an execution query is definitely not syntactically correct, but it gets the point across of what I'm trying to do:

WITH rast AS (
                        SELECT rast FROM "raster2"
                    )
                    CREATE TABLE "raster3" as (select ST_Union(
                        ARRAY[st_mapalgebra(rast, 1, 'callback_fn(float[], integer[], text[])'::regprocedure),
                        st_mapalgebra(rast, 2, 'callback_fn(float[], integer[], text[])'::regprocedure)]
                    ));

It is likely I am way over-complicating the process, but I am fairly new to PostGIS.  One of my concerns is that an approach like this will not preserve the underlying geometry of the raster, which is a must.  How can I change my queries to be functional in the way I need.  Aside from that, am I even approaching this in the right way, or is there a much cleaner way?


texts: How to draw perpendicular/Parallel line using pyQGIS?

texts: While digitizing a line with multiple segment, how to create a particular segment parallel to (or) Perpendicular to another feature using pyQgis?


texts: Filling polygons with patterns in Carto using polygon-pattern-file?

texts:
I'm trying to create a layer in Carto with three categories of transparent hatching of different colors.

I understand how to use the polygon-pattern-file Carto CSS style but am wondering where I can access existing patterns or how to upload my own patterns. It seems like you can't upload .png's to Carto and I can't figure out where else you are supposed to host them in order to get URL's.

I can't seem to find any documentation on this.


texts: QGIS Crashes when selecting > Settings>Project Properties

texts:
I am new to QGIS and this forum but have not found a related problem with a solution.

Running Quantum GIS 1.8.0, on Windows 7 Pro

Running the software, the first action I would normally take on a new project is to set the project properties; however, as soon as I select the project properties option qgis stops working and I get a Windows error reporting notice and QGIS closes.

I have tried re-installing from both the windows installer and the OSGeo4W package and I get the same result every time.

Does Qgis keep an error log?

Any advice would be appreciated.


texts: Difference between iface.mapCanvas and QgsMapCanvas

texts:
I am working on standalone pyqgis script and I wanted to set project extent to position on imported layer.
First I tried with:

vlayer = QgsVectorLayer("data/layers/layer1.gpkg", "layer_name", "ogr")
canvas= iface.mapCanvas()
canvas.setExtent(layer.extent())
canvas.refresh()

And I got: AttributeError: 'NoneType' object has no attribute 'mapCanvas'

Then I tried with using the same layer:

ex = vlayer.extent()
canvas = QgsMapCanvas()
canvas.setExtent(ex)

And then it worked.
I am wondering why i could do this using QgsMapCanvas and not using mapCanvas from iface(which i imported from qgis.utils import iface).

What is the difference between these two object?
Is it meant to be like this because someone would want to work with native Qgis and someone with standalone app or something else?


texts: Configure which bands to display from a 4-band ArcGIS image service?

texts:
The USDA (via the National Map) provides 4-band NAIP ArcGIS "ImageServer" (services) for most states in the US.  I am able to successfully connect to the WMS service in QGIS as described here, which loads the imagery into QGIS using the default RGB (1,2,3) bands.

However, I'd like to change the imagery to display the CIR (bands 4,1,2).  On the "Style" tab, the Render Type only offers the "Singleband color data" option, not the "Multiband color" option described in step 3.5 for typical rasters here.  This appears to be a limitation of the WMS Service.

You can do this from the same service in ArcGIS Desktop using the REST endpoint.  An ArcGIS Online webmap example of changing the bands of a service is available using the REST endpoint.  Just go to the "..." on the layer, choose "Image Display" and then "User Defined Renderer".

UPDATE:

I've been playing around with the Developer Tools in Chrome and the ArcGIS.com map sample from here, and think I've found a small nugget of info.  After I changed the image display settings, I noticed some query parameters set as:

https://gis.apfo.usda.gov/arcgis/rest/services/NAIP/Montana_2015_1m/ImageServer/exportImage?f=image&bandIds=3,0,1&bbox=...

After I changed the band ID's again, I realized that the bandID's in this query parameter are zero-based, instead of 1-based.  Therefore &bandIds=3,0,1 are really referring to bands 4,1,2 used to display imagery in CIR format.  So I just need to figure out how to configure the query string for QGIS to accept this parameter.  Any thoughts?


texts: Reordering different feature classes fom different layers in ArcMap Legend

texts:
I have 6 layers in my TOC.
2 are polygon, 2 are line, and 2 are point feature classes.

TOC looks like this:

Well (water) - point
Pump (Gas) - point
Pipes (Water) - line
Gaspipes (Gas) - line
Gas field (Gas) - polygon
Lake (Water) – polygon

I would like to create Legend in my layout but I need to group different feature classes from different layers.
Legend should look like this:

WATER
X Well
-- Pipes
O Lake

GAS
X Pump
-- GasPipes
O Gas Field

Can this be done without converting legend to graphic?


texts: ArcPy - Force layers to inherit their ColorMap when visualized

texts:
I have a large set of GeoTiff's which after manipulating I apply a colormap to the ones I want to visualize with:

arcpy.AddColormap_management(raster, "", color_map)

I then add them to my MapDocument with:

mapping.AddLayer(df,layer, "BOTTOM")

But they do not automatically use the Color Map Symbology that I had designated but instead use a Stretched Symbology. I want to be able to have the layer use its color map by default either after it is added or before, so that either before I add it or after it's added I can change it to visualize this.

Alternatively, some way to set the imported layer to Apply Unique Values, as this automatically generates the color map manually in the TOC and removes the colors which don't have the values in the color map.

How can this be done? Below is an example of what the symbology should look like "raster01" and what I don't want it to look like (i.e. what it looks like when I import it with AddLayer) "raster02".


texts: Styling points with shapes and graduated symbology in QGIS

texts:
I have a point vector layer based on a .csv file. I am trying to show categorical data (seasons) and continuous data (catch per unit effort, CPUE) at sampling sites that get revisited.

I have used point displacement to account for site revisits, and I have graduated symbology to account for my continuous CPUE data, but now I am having problems styling my categorical data with shapes.

Basic thought outline for styling:
Multiple site revisits   - point displacemet
Continuous CPUE data     - graduated color scales
Categorical season data  - different marker shapes (problem here!).


texts: Can multiple users access ArcGIS styles on server at same time or only one?

texts:
Running ArcGIS 10.0 in our environment.

I created a styles db and saved it on the server. Can multiple users access the styles at the same time or only one?

Second, after experimenting I can add the style with the Style Manager and then click on "Set as Default List" and I think that will make it so that any new map created will reference the styles by default, or is it only for that map?

I don't quite understand how the styles are referenced. Is it per map or per ArcMap session?


texts: Calculating frequency in Field Calculator

texts:
I have a table with many instances of duplicate values.

My end goal is to have three columns:

"The Value": a file name. most are non-unique
"Duplicate Instance": "this is the 1st time this file name appears in the list", the 2nd, the 3rd, etc.
"Total count of Duplicates": "this file name appears 5 times in this list"

I have the  Duplicate Instance from @FelixP's answer to Auto incrementing field based on groups within feature class?

I can create a separate frequency table with the Frequency tool then join that table to calculate in the total count of each file name, but I would love to skip a few steps and directly calculate the frequency in the Field Calculator interface.

Is it possible to directly calculate frequency in ArcGIS Desktop Field Calculator?


texts:
I am importing a gsx set of waypoints with lat/lon taken in degrees, minutes, and thousandths of a minute.  When they are added to my project, the waypoints are not showing up correctly - the coordinates are way off.  New to the software, so appreciate the help.

CH


texts: Workflow for importing a flexviewer website into Flashbuilder 4

texts: I have websites built using the flexviewer, but would like to import them into Flashbuilder 4 in order to enhance them. Is there a way to do this by importing the .xml or config.xml from the compiled flexviewer? Or is the workflow to build custom widgets in flashbuilder and deploy them to the flexviewer?


texts:
I am listing some of the acronyms used in our work and had to stop to think about whether to list these together or separately:

DEM: Digital Elevation Model
DSM: Digital Surface Model
DTM: Digital Terrain Model

In some circles they seem synonymous, in others they seem to have distinct meanings.

Is there an agreed-upon definition for each that enables them to be distinguished?


texts:
I am trying to find bathymetry for Lake Victoria (or portions thereof).

Any GIS format and almost any resolution will do.

As a last resort a hydrographic chart will suffice.


texts: I was wondering if there is a way to make an Autocorrelated DEM within QGIS or Grass.  I am looking to do this with stereo aerial photography.  I am able to do this with a number of COTS tools, but I am looking for an open source solution.  If you know of other tools I am open to suggestions.


texts: How can I output a list of features and # of collected GPS points from Trimble Pathfinder?

texts:
I want to export a list of collected features in an SSF (or .cor) file displaying how many GPS satellite positions were collected for each feature, in order to make sure certain collection standards are being met.  I don't want to select each feature in pathfinder and check the # of positions for each as there are hundreds of features in each SSF file, so a bulk export I could potentially dump into excel to work with would be great.

Is there any way to do this, either as an export from Pathfinder, or as an attribute to add to a feature when they are being exported into a shapefile?


texts:
I am trying to do a temporal analysis on vegetation quality so I downloaded a MODIS .hds raster from earthexplorer.
When I opened the file into ArcMap 10.2.2, it came as a group of 11 rasters.
I stuck on clipping those rasters by my shapefile and combining them into one. (tried, clip raster by shapefile, composite bands...)


texts: Reducing several Points within shapefile to single point

texts:
as an Output of a zonal-statistics calculation I get points representing the lowest points in each polygon. (the Polygons represent Dolines)

How can i reduce the number of Points per Polygon to one?
it doesn't matter which of the Points as they all represent the same value.
In the Attribute-table of the Point feature class the ID of the Polygons is already added.


texts: How create a filter/pattern based on date and time interval between datapoints in QGIS?

texts:
I am working with a dataset of thousands of datapoints in QGIS. The datapoints are satellite pings that are usually taken each 2 hours. But for some reason there are some datapoints that are taken 1 hour or even 15 minutes after each other. I want/need an equally distributed dataset, so I have to filter this data out. There should be a minimum distance of 2 hours between each point.

I have already tried to measure the difference between following datapoints. However, if I would filter on that, I would lose a lot of valuable data because it would filter out any point that is less than 2 hours away from the next point. And there are some parts of the data where hundreds of points just have an interval of 15 minutes. I would lose all of them in this way.

For example:
Point A: 10:15
Point B: 10:45
Point C: 11:15
Point D: 11:45
Point E: 12:15

If I would filter out data based on distance between points <2 hours, I would lose all data. What I want is to create some kind of a pattern that measures =>2 hours ahead from Point A. So it will detect Point E and delete point B, C and D from the dataset.

So does anyone know how to create some kind of a pattern-based filter in the dataset so that there is a minimum time difference of 2 hours between each point? Is anything like that possible in QGIS?

Btw. The data is ordered in this format: YYYY-MM-DD hours:minutes


texts: Querying for points within polygon in MySQL?


texts: Images Shown in Popup Not in Table of ArcGIS Pro

texts:
I'm trying to export data from a map layer to an excel file. From ArcGIS Online I opened the layer using ArcGIS Pro. The images exist in AGOL, and can be seen from the pop-ups in Pro when individual points are clicked. But they aren't in the attribute table for some reason.
Does anyone know how I can get these points to show up in the table as well as the pop-up?


texts: Change in coordinate precision with ogr2ogr moving from GDAL 2 to 3

texts:
I'm migrating some old code from GDAL 2 to GDAL 3. I use ogr2ogr to import a shapefile into PostGIS like this:
ogr2ogr -f PostgreSQL -nln mwe PG:"host=localhost user=myuser" mwe.shp -overwrite

I've saved a copy of an MWE single feature shapefile here if you'd like to follow along. Plain text would make things easier, but the error only seems to occur when I write data from shapefiles to Postgres.
Expected results - GDAL2
When using gdal2, I can query this data in postgres and see the coordinate precision as expected:
select ST_Y(geom)
from (
  select (ST_DumpPoints(wkb_geometry)).geom
  from mwe
) pts;

yields:

Full info for 'working' stack:
ubuntu 18.04
postgres 11
postgis 3
GDAL 2.2.3, released 2017/11/20

POSTGIS="3.1.1 aaf4c79" [EXTENSION] PGSQL="110" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 4.9.3, 15 August 2016" LIBXML="2.9.4" LIBJSON="0.12
.1" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)"

Actually results - GDAL3
When I run this using GDAL 3, some vertices have inflated coordinate precision:

Full info for not working GDAL stack:
ubuntu 20.02
postgres 13
postgis 3
GDAL 3.0.4, released 2020/01/28

 POSTGIS="3.1.1 aaf4c79" [EXTENSION] PGSQL="130" GEOS="3.8.0-CAPI-1.13.1 " PROJ="6.3.1" LIBXML="2.9.10" LIBJSON="0.13.1" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)"

ogrinfo coordinate precision
When running ogrinfo on my MWE shapefile using GDAL 3.0.4, I also see the expected coordinate precision:
ogrinfo -al mwe.shp

INFO: Open of `mwe.shp'
      using driver `ESRI Shapefile' successful.
. . .
 POLYGON ((21.04237 39.579149,21.038807 39.577435,21.036219 39.576187,21.033495 39.574875,21.03046 39.573414,21.024775 39.570475,21.023056 39.569588,21.017654 39.563602,21.009989 39.563313,21.008549 39.563259,21.001377 39.565672,21.001286 39.565702,21.004169 39.575203,21.006231 39.581999,21.013037 39.583546,21.021423 39.584133,21.039852 39.585428,21.04237 39.579149))

Questions

I suspect ogr2ogr is the issue because I'm using PostGIS 3 in both cases - is that reasonable?
I've tried PG_USE_TEXT YES and PRECISION=NO flags in ogr2ogr to no effect
Is there any chance that the higher precision coordinates are what's stored in the shapefile, even though ogrinfo with GDAL 3 doesn't report them?

I could certainly truncate / manage this change in PostGIS, but am trying to understand why this happened and if there's a way to get the previous behavior.


texts:
I am attempting to identify utility poles that require inspection. Currently, the feature has a field containing the pole's birthmark in a string format and an installation date field (with the majority being null values at the moment). The birthmark field is formatted as 40-5-18 where 40 is the height, 5 is the class, and 18 is the year. I need a Python solution to capture the last two digits of the birth mark string and either update the installation date field to some arbitrary date in the corresponding year or select values (< or > a particular year) using a mathematical operator. The ultimate goal is to be able to query poles due for inspection based on the date of the birthmark or installation.

Thoughts?


texts: Issue with address locator in ArcGIS 10.1

texts:
I am attempting to geocode a table of addresses for a client. I researched into how to create your own address locator and have two separate reference sources: one that has exact address information for our county and another that has road information.

I have made a locator for addresses using the US Address - Single House style, and one for roads (I experimented with one range/dual ranges style). When I make a composite locator using both of these and add it to ArcMap during geocoding, I am not presented with an "Address" field input. I only have street or intersection, city, state, zip.

What setting exactly will allow me to input the physical house number?

Edit:
I should mention that the client data has addresses split into multiple fields such as: House_NO, ST_NAME, CITY.

Should I simply combine the house_no and st_name fields?


texts: Mask clouds across multiple images using Google Earth Engine

texts:
I Have multiple Landsat8 images at the same place in different years.

I want to have a cloud mask for each image, then merge the mask so that all images will be masked with a same aggregated mask.

I've got masked image for each year using code from this post Using cloud confidence to create cloud mask from Landsat 8 BQA?
.

var RADIX = 2;  // Radix for binary (base 2) data.
var extractQABits = function (qaBand, bitStart, bitEnd) {
  var numBits = bitEnd - bitStart + 1;
  var qaBits = qaBand.rightShift(bitStart).mod(Math.pow(RADIX, numBits));
  return qaBits;
};

var getMaskedImg = function (i){

  // Reference a sample Landsat 8 TOA image.
  var image = ee.Image(i).clip(geometry);

  // Extract the QA band.
  var image_qa = image.select('BQA');

  // Create a mask for the dual QA bit "Cloud Confidence".
  var bitStartCloudConfidence = 5;
  var bitEndCloudConfidence = 6;
  var qaBitsCloudConfidence = extractQABits(image_qa, bitStartCloudConfidence, bitEndCloudConfidence);
  // Test for clouds, based on the Cloud Confidence value.
  var testCloudConfidence = qaBitsCloudConfidence.gte(2);

  // Create a mask for the dual QA bit "Cloud Shadow Confidence".
  var bitStartShadowConfidence = 7;
  var bitEndShadowConfidence = 8;
  var qaBitsShadowConfidence = extractQABits(image_qa, bitStartShadowConfidence, bitEndShadowConfidence);
  // Test for shadows, based on the Cloud Shadow Confidence value.
  var testShadowConfidence = qaBitsShadowConfidence.gte(2);

  // Calculate a composite mask and apply it to the image.
  var maskComposite = (testCloudConfidence.or(testShadowConfidence)).not();
  var imageMasked = image.updateMask(maskComposite);
  return(imageMasked)
}

How can I merge the masks and let all images share this merged mask?


texts: Is there any map API that allows you to pull an address range from a polygon?

texts:
I've seen a lot of answers saying you can't draw polygons and get the exact addresses of all of the places contained; thankfully, I don't need so much info.

I'm looking for a way to get a vague address range from within a set of coordinates. Say, draw a triangle and see that it encompasses..

2nd street, 120 - 280;

3rd street, 80 - 320;

4th street, 40 - 460;
 and so on.

My goal (since there's probably a much simpler solution I'm not thinking of) is to be able to have a database of address ranges attributed to, say, municipalities in Excel. When typing an address in, a row could auto-populate applicable municipality name (or country, postal/zip code block, local sports league, school jurisdiction, whatever shapes you've defined) based on the range info.

A1"110" B1"3rd" C1"STREET" is within polygon: [Sheet2A1]"crazy ex, avoid" range at [Sheet2Betc.]'3rd street' <80, 320>.

D1 Populates with [Sheet2A1]


texts: Calculate shortest path not on land but over sea

texts: I wonder if someone can give me some advice on how to calculate the shortest path between two points, not over land but over sea. I know that I will have to take in to account the magnitude of sea currents for specific periods of time, but besides that I am a little bit lost.


texts: Rings not displaying on .TAB import to .SHP save

texts:
I have some strange behaviour happening when importing a MapInfo .TAB file into QGIS. One feature contains internal rings which are ok on initial import but do not show on the first .SHP save. I can see the ring perimeter but selecting picks the feature as a whole, shading in the ring areas. It's only if I save again that they appear properly cutout and all appears ok.

I've used the MapInfo file with Manifold and no issues in that software. However saving the file as a .MIF in Manifold and importing into QGIS I lose the rings in the initial import and again need to save twice to .SHP to get them to appear ok.

Just wondered if anyone had an explanation?


texts: Arcpy Subtracting rasters error - TypeError: 'unicode' object is not callable

texts:
I have a script which subtracts 2 versions of the same raster from each other (one newer than the other) and then extracts the additions & deletions that have happened.

I created a simple test script which works:

# Import arcpy module
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Local variables:
New = "C:\\Data\\JOB_13\\tq03\\ABRdsText.tif"
Old = "D:\\Data\\Job_Backups\\JOB_13\\tq03\\ABRdsText.tif"
Addition_Raster = "C:\\Testing\\District_Editor\\QC_ChangePlots\\JOB_13\\tq03\\Additions\\output.tif"
Deletion_Raster = "C:\\Testing\\District_Editor\\QC_ChangePlots\\JOB_13\\tq03\\Deletions\\output.tif"

arcpy.MakeRasterLayer_management(New, "N")
arcpy.MakeRasterLayer_management(Old, "O")

# Process: Raster Calculator
RR = Raster("N") - Raster("O")

# Process: Extract by Attributes
arcpy.gp.ExtractByAttributes_sa(RR, "VALUE =-255 OR VALUE =0", Addition_Raster)

# Process: Extract by Attributes (2)
arcpy.gp.ExtractByAttributes_sa(RR, "VALUE =255 or VALUE=0", Deletion_Raster

I then replicated this into my main script:

#import the necessary modules
import arcpy
from arcpy import env
from arcpy.sa import *
import os

#Checkout the necessary Arc licences
arcpy.CheckOutExtension("spatial")

#Ask the user for the job number
#Job_Number = arcpy.GetParameterAsText(0)
# MANUAL OVERIDE TO GET TEST JOB_NUMBER
Job_Number = 13
Job = 'JOB_'+str(Job_Number)

####################################################
#Folder where all data will be written to
OutputsDirectory = r"C:\Testing\District_Editor\QC_ChangePlots"

####################################################
#Folder with the backed up data pre-editing i.e. old data
BackupDirectory =r"D:\Data\Job_Backups"

####################################################
#Location of the QC Change plot Master mxd
ChangePlot_mxd = r"Y:\esriuk\OSBurnProcessGPTools\resources\mxdHolder\QC_ChangePlot.mxd"

####################################################
#Location of the Job folder
DataFolder = os.path.join("C:\\Data\\",Job)

####################################################

JobF = os.path.join(OutputsDirectory,Job)

for folder in os.listdir(DataFolder):
    arcpy.env.workspace = DataFolder + "\\" + str(folder)
    Ras = arcpy.ListRasters("*", "TIF")
    for Raster in Ras:
        #Set the output rasters
        Additions = os.path.join(JobF + "\\" + str(folder) + "\\Additions\\"+ str(Raster))
        Deletions = os.path.join (JobF + "\\" + str(folder) + "\\Deletions\\"+ str(Raster))

        #Set the file paths for the new and old rasters
        New = os.path.join(str(DataFolder) +"\\" + str(folder) + "\\" + str(Raster))
        Old = os.path.join(str(BackupDirectory) + "\\" + str(Job) + "\\" + str(folder) + "\\" + str(Raster))

        #Convert the Raster datasets into a temporary Raster Layer so that they can be used for a calculaion
        arcpy.MakeRasterLayer_management(New, "N")
        arcpy.MakeRasterLayer_management(Old, "O")

        #Raster Calculator - subtract the old raster from the new raster ( new - old)
        CalculatedRast = Raster("N") - Raster("O")

        #Extract the Additions
        arcpy.gp.ExtractByAttributes_sa(CalculatedRast, "VALUE =-255 OR VALUE =0", Additions)

        #Extract the deletions
        arcpy.gp.ExtractByAttributes_sa(CalculatedRast, "VALUE =255 or VALUE=0", Deletions)

However, I get the error

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.1\Lib\site-packages\Pythonwin\pywin\framework\scriptutils.py", line 325, in RunScript
    exec codeObject in __main__.__dict__
  File "S:\Development\Dan\GIS\Code\District_Editor\QC_ChangePlots.py", line 204, in <module>
    CalculatedRast = Raster("N") - Raster("O")
TypeError: 'unicode' object is not callable

Anybody able to suggest where I am going wrong and how to fix it?  I'm fairly new to python, so it's probably something quite simple!

The error seems to be here

CalculatedRast = Raster("N") - Raster("O")


texts:
I'm loading a WMS map using gdal.

Online access works fine, including cache. The cache looks fine, I can see the files and open the images.

When I try to load the same WMS in offline mode it works, but only if I use the same coordinates.

My issue is to open the map in offline mode for different coordinates that overlaps with the cached area. (It's ok if the area outside my cache would be black)

Sample Xml for the request as described in the gdal doc, the problem can be reproduced changing OfflineMode to true and changing the coordinates:

<GDAL_WMS>
    <Service name="WMS">
        <Version>1.1.1</Version>
        <ServerUrl>https://ows.terrestris.de/osm/service?</ServerUrl>
        <Layers>OSM-WMS</Layers>
    </Service>
    <DataWindow>
        <SizeX>500</SizeX>
        <SizeY>500</SizeY>
        <UpperLeftX>9.892920</UpperLeftX>
        <UpperLeftY>50.626900</UpperLeftY>
        <LowerRightX>13.978200</LowerRightX>
        <LowerRightY>47.082800</LowerRightY>
    </DataWindow>
    <Cache>
        <Path>./gdalwmscache</Path>
        <Depth>2</Depth>
        <Extension>.jpg</Extension>
    </Cache>
    <OfflineMode>false</OfflineMode>
</GDAL_WMS>


texts: ArcGIS Pro Field Calc to calculate pipe volume

texts:
I'm trying to do this calculation:

(!DIAMETER!/2)*( !DIAMETER!/2)math.pi/144 !Shape.STLength()!*7.48

and keep getting an error.

Can anyone give me the correct syntax for calculating the volume of a water pipe?


texts: Atlas map generation in QGIS using point layer only

texts:
I am trying to automate generating maps with Atlas using points only. As the snapshot shows the data points for each village is shown by a color. Every four villages fall under one health area. There are more than 100 health areas in my dataset. I am trying to create one map per health area. Is that possible?
I know that there is one way doing this by constructing a polygon for each health area, but I rather not to do it because there is no specific border between some villages and they overlap in some areas.


texts: ArcGIS Web Layer a snapshot or real-time proxy

texts:
I'm trying to understand if a Layer generate from the web in ArcGIS is a static snapshot, a proxy for the underlying API or maybe polls.

From here I can add from a WFS API or a CSV file.  I could probably do either, but I do need the ArgGIS layer to stay in sync with changes I make to the underlying data.

Is this how it works, and if not, is it possible to proxy from something other than another Esri source?


texts: Distance to coastline in meters using Cartopy and Shapely

texts:
I am trying to calculate the distance to coastline for some points. Using Cartopy and Shapely, I can calculate the distance in degrees, but that is not very helpful. Perhaps if I could transform the geometries, this would not be difficult, but I'm not sure how to do that in Cartopy. I'm trying to wean myself off the depreciated Basemap package. Does anyone know how to do this where the distance is output in meters? I have a dataframe with lat and lon for the points I want to calculate.
import numpy as np
import pandas as pd
import shapely
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

land = shapereader.gshhs(scale='h', level=1)

geoms = list(itertools.chain.from_iterable(geom.geoms for geom in shapereader.Reader(land).geometries()))
geometries = shapely.geometry.MultiPolygon(geoms)

df['dist'] = np.nan

for i in df.index:
    point = shapely.geometry.Point(df.iloc[i]['lon'], df.iloc[i]['lat'])
    df['dist'][i] = geometries.distance(point)

print(df)


texts: Transferring calculated results from spatial join to original feature class in ArcMap

texts:
Is there a way to transfer the attributes calculated in a field within a spatial join back into the original feature class using ArcMap?

I am numbering manholes ID's based off a grid system, so each grid cell would have manholes with ID's from 1-n. Each manhole ID also has the grid cell name attached to each ID. Ex:

Grid cell J10:
 Manholes: J10-1, J10-2, J10-3, etc...
Grid cell A5:
 Manholes: A5-1, A5-2, A5-3, etc...

The new ID's are stored in a field called New_ID. After calculating those new ID's, I noticed that the original New_ID field in my original feature class had not updated like it would with a tabular join. And I can't copy and paste my results from the spatial join field into my original feature class New_ID field. Does anyone know if there is a way to update a field in your original feature class using data you calculated in a spatial join?

I am using ArcMap 10.4.


texts: Does anyone know what a geoid10 and its numbers mean?

texts: I am doing home for my GIS class where we had to download census data from our county and to define what geoid10 and its numbers mean. I cannot find the definition anywhere.


texts: Buffer zone outside objects cluster using FME

texts:
I am doing a buffer analysis around some roads. The buffer is 1000 meters and is dissolved. This gives me a buffer which covers the roads as well as the distance outwards. It looks like this:

Instead of having a full buffer also covering the area in between the roads I would like to have a 1000 buffer outside the roads.

I was thinking I could do a buffer along the outer line of the buffer of 1000 meter and then do an AreaonAreaOverlayer where I select the polygons which overlap. Could work - but if there is an other way using FME I would be thrilled.


texts: ArcPy - Select Layer by Attribute / Delete Selected Layer - SQL

texts:
I feel like I've tried so many iterations of this, can anyone assist with how to use ArcPy to selected layer by attribute then delete selected layer? I think my SQL query is incorrect?

 import arcpy
 import pandas as pd
 import numpy as np

 pid_field = ["pid"]
 target = arcpy.da.TableToNumPyArray(Master_Subjects, pid_field)

    Master_Subjects = r"C:\Users\dsser\Documents\ArcGIS\Projects\MyProject1\Master_Subjects.shp"
    target_pid = target[0]["pid"]
    # just selecting a random pid here from the Master_Subjects shapefile to see if I can get the select by attribute to work.This is a shapefile feature class with point geometry.
    # Also, target_pid is an integer
    target_pid = target[0]["pid"]
    query = "'" "subject_pid" = target_pid "'"
    arcpy.SelectLayerByAttribute_management(Master_Subjects, "NEW_SELECTION", query)
    arcpy.DeleteFeatures_management(Master_Subjects)


texts: Defining map projection of Eastern part of Ukraine

texts:
I have a map of the Eastern part of Ukraine, which was apparently used in the ceasefire negotiations to define the borderline between the conflicting sides.

The full map: https://informnapalm.org/wp-content/uploads/2015/01/karta-polnayacc.jpg

The map is quite detailed but I could not find the type of projection it uses.
Moreover, it is difficult (though possible) to see the exact numbers.

However, there are some numbers.

Coordinate X varies from 7288000 to 7612000
Coordinate Y varies from 6192000 to 6608000
The gap between meridians and parallels is equal to 4000 map units.

Coordinates X and Y should roughly correspond to

Longitude 47.12°N - 49.65°N
Latitude 36.15°E - 40.19°E

in WGS84 respectively (I used Google Maps).

The middle of the map is a city called Debaltseve located at 48.33°N and 38.36°E.

Although the headers are in Russian, the names of the places are in Ukrainian.

Thus, I thought that it is likely to be one of the Ukrainian projection systems. I tried to use the epsg.io website and play with the common Ukrainian projection systems it stores but every time I tried to transform the WGS84 coordinates into one of the Ukrainian projection systems, I got an error.

Could you suggest me what kind of projection system the map might use?


texts: Connecting to ArcSDE with ArcSDE Geoserver extension?

texts:
I have a server running Geoserver 2.8 with the ArcSDE extension installed. I have added a new "ArcSDE ESRI(tm) ArcSDE 9.2+ vector data store" with credentials from our DBA.

When I try to add a layer I only see a tiny subset of the layers that should be there. Talking with the DBA it appears the majority of the layers are prefixed with "SDE_" and it is these layers that are not showing up.

We are not connecting with the SYSTEM, DBA or SDE users. I know the layers are prefixed with "SDE_" because when talking to the database guy here he said "We just name them like that so we know they're spatially enabled rather than just tabular when looking through the schema listing"

Does anybody know why these layers are not able to be loaded by the extension? Is there a way of making the extension load them?

Database: Oracle Standard Edition One
Version: 11.2.0.2
ArcSDE Version: 10.0
ArcGIS Geodatabase Version: 10.0 for Oracle11g Build 2004
SDE Client library: jpe_sdk.jar & jsde_sdk.jar (Created-By: 1.4.2_12 (Sun Microsystems Inc.), Specification-Version: 9.3)
JDBC Driver version: ojdbc6.jar (Implementation-Version: 11.2.0.4.0, Repository-Id: JAVAVM_11.2.0.4.0_LINUX.X64_RELEASE)
Apache Tomcat version: 7.0.65


texts: Converting wkt to geojson with OpenLayers


texts:
I asked a question earlier about how to invert my file "hydropolys.shp" from a vector of ocean and seas to a vector of land and got sent to here, but I am experiencing issues where using the "Difference" tool, it returns an empty map.

So I tried to follow the instructions geoprocessing but this tool keeps giving the "Not Responding" error. Any suggestions on what to do next? Or any other easy way to convert the shapefile? I am using QGIS 2.18. Or is there a way to do this instead in code in PyQGIS?


texts:
I think I understood the steps to write an array to a raster file using python. however it doesn't work in my case and I don't know why. I have an array that I want to write to a Geotiff with Swiss Oblique Mercator Projection. I have a DEM as Geotiff in the exact same location with same projections and Geotransform, so I thought I just read all parameters from that DEM and use them to write my new raster.
Here's the code that I'm, using:

driver.Register()
fn = '/media/LACIE_SHARE/davos2015.tif'
inDs = gdal.Open(fn)
if inDs is None:
    print 'Could not open ' + fn
    sys.exit(1)

rows=inDs.RasterYSize
cols=inDs.RasterXSize

driver=inDs.GetDriver()
outDs=driver.Create('/home/grassdata/test.tif',cols,rows,1,gdal.GDT_Float32)

if outDs is None:
    print "crap"
    sys.exit(1)
outBand=outDs.GetRasterBand(1)
outData=coh0 #coh0 is the array I want to convert to raster. it's an array with float32 entries.

outBand.WriteArray(outData) # this gives 0
outBand.FlushCache()
outBand.SetNoDataValue(-99) #this gives 0
outDs.SetGeoTransform(inDs.GetGeoTransform()) #this gives 0
outDs.SetProjection(inDs.GetProjection()) #this gives 0
del outData

the lines where I commented accordingly raise 0 as result and I don't understand why.


texts: How to get the wavelength values of bands from hyperspectral image?

texts:
I am trying to open and read a hyperspectral image using GDAL like below.

So after reading the data like this :

data = gdal.Open( filename )

then I can see the x, y, dim information of the dataset :

x = data.RasterXSize
y = data.RasterYSize
dim = data.RasterCount

Now the question is that if there is any way to find out that each band wavelength from image.?
I have the following data image:
sub_66 Type:File
sub_66 Type:hdr
sub_66 Type:aux(xml)
This data file is of AVIRIS(Hyperspectral Data) sensor.


texts: How to create multiring buffer with sf in R

texts:
I have a polygon of an area which I like to buffer multiple times with different distances for intersection with another layer buffer(ring)-wise later (this part already works!)

What I do not achieve, is to create a single sf object from the list of buffers (to be able to export it as a layer).
The intersections from another  layer with those buffers (done in the same loop) DOES get combined to a sf object by do.call(rbind ,...)!

I'm lost in the complexity of those object types.

Example:

#Polygon
x <- structure(list(AREA = 1489841.35074, geometry = structure(list(
    structure(list(structure(c(4371395.39166986, 4371395.39166986,
    4372634.07252284, 4372634.07252284, 4371395.39166986, 5478810.68590502,
    5480013.45040013, 5480013.45040013, 5478810.68590502, 5478810.68590502
    ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"))), n_empty = 0L, crs = structure(list(
    epsg = NA_integer_, proj4string = "+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +datum=potsdam +units=m +no_defs"), .Names = c("epsg",
"proj4string"), class = "crs"), class = c("sfc_POLYGON", "sfc"
), precision = 0, bbox = structure(c(4371395.39166986, 5478810.68590502,
4372634.07252284, 5480013.45040013), .Names = c("xmin", "ymin",
"xmax", "ymax"), class = "bbox"))), .Names = c("AREA", "geometry"
), row.names = c(NA, -1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(NA_integer_, class = "factor", .Label = c("constant",
"aggregate", "identity"), .Names = "AREA"))

library(sf)

#distances for buffers                                                                                "aggregate", "identity"), .Names = "AREA"))                                                                                "aggregate", "identity"), .Names = "AREA")
dist <- c(50,500,1000,3000)

#list holding the buffer rings
lstbuffers <-list()
for (i in seq(dist)){
  #first buffer
  if (i == 1){
    lstbuffers[[i]]<- st_buffer(x,dist[i])

  }else{
  #next buffers as rings around the previous
    lstbuffers[[i]] <- st_difference(st_buffer(x,dist[i]),st_buffer(x, dist[i-1]))
  }
  lstbuffers[[i]]$radius <-dist[i] #add the radius
}

buffers <- do.call(rbind,lstbuffers) #DOES NOT WORK, I get a matrix of lists :(

So, I can do operations with my lstbuffers[i] list items, but fail to combine them into one layer.


texts: Where can I find annual temperature data for every city and town in the US and Australia?

texts:
Where can I find yearly annual averaged temperature data for every city and town in the US and Australia for a heatmap?

links are good


texts:
Why are otherwise valid geography polygons (st_isvalid(col::geometry) is true) invalid when trying to calculate the area, when they lie on the equator?

For example:

select st_area(st_geographyfromtext('POLYGON((-148.282525588539 0,-147.881585011895 0.29129549292729,-148.43567339643 0.471329293759851,-148.282525588539 0))'))

The above returns

ERROR: ptarray_area_spheroid: cannot handle ptarray that crosses equator


texts: Converting ENVI file to TIFF using ArcPy

texts:
I have ENVI files one ASCII file .hdr and a flat_binary raster file which is a file without an extension
I'd like to convert them to TIFF files
These are the steps I followed:

Rename flat_binary raster file by adding “.dat”
Keep the pair of files in the same directory
Then, I wrote this code

import arcpy
ndvi_file="G:\\data\\test\\MYD13A2_A2006201_MED_1km.dat"
ndvi_tiff_file="G:\\data\\test\\data_tiff"
arcpy.RasterToOtherFormat_conversion(ndvi_file, ndvi_tiff_file, "TIFF")

But I didn't get the converted File
I'm getting this error
Failed to convert: <built-in function input>. ERROR 999999: Error executing function.
Failed to execute (CopyRaster)

Is there anything missing?


texts: How to change default annotation properties in ArcMap 10.4?

texts:
I am trying to figure out how to change the default properties of an annotation feature class in ArcMap 10.4. I see the "Text Symbol" box when I double-click on the "Default" category under my annotation feature class in the table of contents, but every option in it is greyed out.

Also, I cannot figure out how to make a new template for annotation either. I know how to do it with regular features in the Create Features window. Any advice would be welcome.


texts:
The recent release of "new Google Earth" appears completely unsuited to the tasks for which I use Google Earth, so I figured it was time to see if I should migrate to ArcGIS Earth as was suggested to me last year.

I installed ArcGIS Earth and found a way to import My Places.kmz from Google Earth.  However it acts like all the data is read-only.  I cannot figure out how to edit placemarks and paths or even rearrange them in the Layer List.

Is it possible to edit imported KML in ArcGIS Earth?


texts:
Does my code have problem? Because if I run in Visual Studio, it doesn't show error, but the file didn't load on my map.
I put the JSON file and HTML on my desktop.

This is my code :

    var osmUrl = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
        osmAttrib = '&copy; OpenStreetMap contributors',
        osm = L.tileLayer(osmUrl, { maxZoom: 20, attribution: osmAttrib });
    var map = L.map('map').setView([24.151687694799833, 120.64116954803465], 15).addLayer(osm);

    L.control.scale().addTo(map);

    $.getJSON("taiwan.json", function (data) {
        var datalayer = L.geoJson(data, {
            onEachFeature: function (feature, featureLayer) {
                featureLayer.bindPopup(feature.properties.date);
            }
        }).addTo(map);
        map.fitBounds(datalayer.getBounds());
    });

</script>

This is my json (part of one)

{"type": "FeatureCollection","features": [
{
  "type": "Feature",
  "properties": {
    "date": "1\/5\/1995 6:14:00",
    "id": "第001號",
    "lat": "24.96",
    "lng": "121.7",
    "depth": "91.5",
    "magnitude": "4.5",
    "desc": "台北市地震站東偏南方20.2公里",
    "_id": 49710583
  },
  "geometry": {
    "type": "Point",
    "coordinates": [ 121.7, 24.96 ]
  }
},
{
  "type": "Feature",
  "properties": {
    "date": "1\/10\/1995 15:55:00",
    "id": "第002號",
    "lat": "23.68",
    "lng": "121.43",
    "depth": "3.8",
    "magnitude": "5.1",
    "desc": "花蓮西林地震站南方14.6公里",
    "_id": 49710584
  },
  "geometry": {
    "type": "Point",
    "coordinates": [ 121.43, 23.68 ]
  }
},


texts:
I'm attempting to create a series of maps in QGIS atlas with the same coverage and background data but different raster outputs (flood mapping outputs for a range of scenarios).

I've created a shapefile coverage layer containing polygons for each of the scenarios and using the polygon identifier attempted to turn on and off raster layers using a function within a label object. I've played around with the function for a while but can't get it to work...

Does this sound like the right approach that I'm taking? Is there a better way to do this?

I'm very much a beginner when it comes to using functions in QGIS.


texts:
When using a cursor to loop through an attribute table, is it possible to iteratively select each row?

I want to clip a line to a polygon grid (see previous question); due to the size of the grid the split function fails. Hence I'm interested in selecting each polygon, and clipping the line to the selected cell by using the shapefile with only that cell selected.

I'm familiar with looping through the attribute table and performing calculations - but I don't know how to select a row.


texts:
I am trying to do the following calculation:

Given lat,lon of a point (center of a rectangular (in UTM) domain), I find its UTM coordinates.
Given domain width and length, I find coordinate of its 4 corners in UTM
Finally, I need to find lat,lon coordinates for the corners

I'm having a problem with step 2, because one of the corners is outside of the UTM zone in which the center point is located. So I get negative easting numbers and can't figure out how to do the lat,lon conversion. I'm supposing that I need to find the new UTM zone for that point, but I don't know exactly how. Any help will be appreciated.

Additional info: I'm using this little Python library to do the UTM<->lat,lon conversion https://github.com/Turbo87/utm It doesn't seem to handle negative eastings.


texts: Simultaneous multiple vector layer intersection algorithm

texts:
Is there an algorithm, you can suggest or recommend, that I can implement to intersect multiple vector layers?

I have a few hundred layers I want to intersect and thinking of developing a little QGIS python script to solve it (it doesn't have to be QGIS). All the layers are of the same geometric type of either polyline or polygon.

I am aware that I can intersect them two at a time but perhaps it's more efficient to implement a mechanism that intersects all these layers at the same time?


texts:
I am using geomintersects in QGIS and there is, unfortunately, a limit for feature test of this function (100000  test).

Is there a way to exceed this limit to take more than 100000 test?


texts: How to rasterize an object of class SpatialPolygonsDataframe in R

texts:
When I read a polygon shapefile into R using rgdal, it was imported as an object of class SpatialPolygonsDataframe.

I need to rasterize the polygon, but the rasterize{raster} function requires it to be of the SpatialPolygons class.

How do I convert it from the current data frame into the SpatialPolygons class?


texts: How delete feature from layer using Gdal binding java

texts:
I am trying delete one feature from layer using GDAL in java https://trac.osgeo.org/gdal/wiki/GdalOgrInJava

But nothing happens.

    ogr.RegisterAll();
    Driver driver = ogr.GetDriverByName("ESRI Shapefile");
    DataSource ds = driver.Open(shape,1);

    Layer layer = ds.GetLayer(0);
    System.out.println(layer.GetFeatureCount());

    layer.DeleteFeature(0);
    ds.SyncToDisk();

Anyone know what the problem is?
I tried too use layer.StartTransaction(); and layer.CommitTransaction();


texts: PostgreSQL: transform spatial column from PG_GEOMETRY to ST_GEOMETRY

texts:
I have spatial Data stored in my PostgreSQL database as PostGIS Geometry (PG_GEOMETRY). It looks for example like this: 010100002031BF0D0000000000F989264152B81E057E6D5841

Now I would like to move this column to a database with ArcSDE schema, which uses ST_GEOMETRY.

Is this possible with a simple sql command, or a more difficult task?


texts: Does WMS GetFeatureInfo takes into account the applied Style (SLD)?

texts:
As the title says, does WMS GetFeatureInfo takes into account the applied style to a layer? For example I have layer of all US States and I applied the following SLD:

<PropertyIsEqualTo>
   <PropertyName>STATE</PropertyName>
   <Literal>California</Literal>
</PropertyIsEqualTo>

so that only the State of California will be visible. When I click randomly outside California, should it still give the info of the other states? If so, is that the default behaviour? Is there any way of controlling this?

Thanks a lot!


texts:
I've just changed to QGIS 3 and I can't figure out why I can't change from Python 2.7 to Python 3.6 in PyCharm, keep getting invalid SDK!
I run this .bat-file:

SET OSGEO4W_ROOT=C:\OSGeo4W64
SET QGISNAME=qgis
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
SET QGIS_PREFIX_PATH=%QGIS%
SET PYCHARM="C:\Program Files\JetBrains\PyCharm 2017.1.2\bin\pycharm64.exe"

CALL %OSGEO4W_ROOT%\bin\o4w_env.bat

SET PATH=%PATH%;%QGIS%\bin
SET PYTHONPATH=%QGIS%\python;%PYTHONPATH%

start "PyCharm aware of QGIS" /B %PYCHARM% %*


texts:
import arcpy
rows = arcpy.SearchCursor("E:/Arcgis/lx.gdb/testSUBorigin", "", "", "",  "VEHICLEID A; TIME A")
arcpy.AddField_management("testSUBorigin", "SUB", "TEXT")
cursor = arcpy.SearchCursor("E:/Arcgis/lx.gdb/testSUBgetondropoffshp")
for line in cursor:
    for row in rows:
        if (line.VEHICLEID == row.VEHICLEID) and (line.TIME != row.TIME):
           row.SUB = "36_1"
         else:
           line.next()

  Runtime error  Traceback (most recent call last):   File "", line 3, in    File "d:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 3243, in AddField     raise e ExecuteError: ERROR 000732: Input Table: Dataset testSUBorigin does not exist or is not supported

I want to add Subtrack point name attribute according to the get_on and drop_off points,But I don't know the mistake,and I don't know whether the idea is true.I sincerely hope someone do the similar work can help me.