texts: Use QgsAbstractProviderConnection in PyQGIS (3.14)

texts:
I created a PostGIS connection named 'test_conn' and try to access it via
from qgis.core import QgsAbstractProviderConnection
con = QgsAbstractProviderConnection('test_conn')

as described in https://qgis.org/pyqgis/master/core/QgsAbstractProviderConnection.html?highlight=qgsabstractproviderconnection#module-QgsAbstractProviderConnection
but I get
Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
TypeError: qgis._core.QgsAbstractProviderConnection represents a C++ abstract class and cannot be instantiated

Is this intended behavior or a bug or am I missing something?


texts: Converting spatial data to table of vertices using PostGIS?

texts:
I have a spatial data table in PostGIS 2.1 that looks like this:

Geom                | type |  cost
POLYGON(0 0,1 1,2 2)|  A   | 123
POLYGON(3 3,4 4,5 5)|  b   | 123

I am trying to create a view of this table that looks like this for a system that does not support spatial data natively:

vertexorder  |  x  |  y  | type  |  Cost
1            | 0   | 0   |   A   | 123
2            | 1   | 1   |   A   | 123
3            | 2   | 2   |   A   | 123
1            | 3   | 3   |   b   | 123
2            | 4   | 4   |   b   | 123
3            | 5   | 5   |   b   | 123

Is there anyway to query for vertex like this?


texts: How to create boundary polygons of the actual data area in LAS files?

texts:
Using command line utility or utilities suitable for scripting how can I generate a boundary polygon that shows the area actually covered by lidar data points, not just the box extent of the LAS files?

In this figure, the purple is a simple box extent, easily generated by several tools, but what's desired is the grey area which represents the actual point coverage. Using pdal tindex without --fast_boundary generates a useful but fairly crude approximation. I'd like something that matches the data area more tightly (a bit inside even better).

.

I'm using Windows but don't want to restrict solutions to a single platform unecessarily. Open source is preferred, then free, then whatever.

Related question for ArcGIS:
Creating boundary polygon shapefile from set of LAS files using ArcGIS Desktop?


texts: How to insert "Topo North Arrow" using ArcObject

texts: Does anyone know how to insert "Topo NorthArrow" using ArcObject C# SDK? There is no example or snippet. I am not sure if it is even available through SDK.


texts: SQL Error Message: Find and Replace for Date Field

texts:
I am attempting to use the find and replace function within the Attribute Table Options to replace fields that contain dates in the following format: month/day/year (e.g. 03/14/2014). The field's properties type is DATE.

However, when I try to do this, I receive the following error message: Search Error, an invalid SQL statement was used.

I realize you can set up a field calculator expression to find and replace DATE strings but is it possible to do this within the Attribute Table Options Find and Replace or does the special character “/” in the DATE field preclude this?
See image for ArcGIS  version.


texts: How to Auto Increment Attribute Based on Nearest Point?

texts:
An authority I am working with would like their Hydrants, a point shapefile, to be renumbered based on their location. I am looking for any python script that would help me A. Find the next closest Hydrant B. Generate a number C. Increment that number for all hydrants.

I am using ArcGIS Desktop version 10.0.


texts: Finding points within various radius's of another point

texts:
I am trying to find the number of crash accidents locations that lies within the radius's 0m-100m, 100m-200m, and 300m-400m of each red light camera location. However, I am not sure how I can do this.

I have read a number of ways online on how to find points that lie within a specific radius of another set of points, but I'm not sure how to find points that lie within an area radius around the point. For example, the number of accidents between 100 meters to 200 meters of a red light camera location. This would not include the number of accidents between 0-100 meters of the camera.


texts: Where can I find Tampa Bay building footprint GIS data?

texts: I have looked on the city and county website and can't seem to find it.


texts: Create additional 'fields' in GDB domain system table (GDB_items)?

texts:
Is it possible to create additional 'fields' in the geodatabase domain system table (GDB_items) in SDE 10.3/Oracle?

It would be handy to store additional information in domains.
For example:

CODE [domain coded value]
NAME [domain description]
Additional fields: DETAILS, COMMENTS, RELATED_DOCUMENT, etc.

I know it wouldn't be accessible in ArcGIS desktop's out-of-box functionality, but maybe it would be possible to query it with SQL?

The SQL used to access the domain would look like this, but include additional 'fields':

SELECT
  EXTRACTVALUE(CodedValues.COLUMN_VALUE, 'CodedValue/Code') AS DOMAIN_CODE,
  EXTRACTVALUE(CodedValues.COLUMN_VALUE, 'CodedValue/Name') AS DOMAIN_DESCRIPTION
FROM
  SDE.GDB_ITEMS_VW
  JOIN SDE.GDB_ITEMTYPES
  ON SDE.GDB_ITEMS_VW.Type = SDE.GDB_ITEMTYPES.UUID,
  TABLE(XMLSEQUENCE(XMLType(Definition).Extract('/GPCodedValueDomain2/CodedValues/CodedValue'))) CodedValues
WHERE
  SDE.GDB_ITEMS_VW.NAME = 'STRATEGY'

Or, as an alternative, would it be an acceptable practice to create a separate table that is related to the domain by the CODE value? And then join to it using the query mentioned above?


texts: Getting same results from Mapping toolbox of MATLAB as ArcMap?

texts:
I am newbie in GEO. I have a point [lat,lon] = [-77.5100 , 166.4000]. I want to find its projection based on equidistant Azimuthal projection method. When I use ArcMap 10.3.1 (with default parameters), I receive [x,y]=[1877942.02210000 , 3755681.82180000]. however, when I use MATLAB (origin set to [0,0]), I receive [x,y]=[ 0.0927, -1.7802]. How can I tune MATLAB to produce the same result of ArcMap.

Details:
MATLAB code:

mstruct = defaultm('eqdazim');
mstruct.origin = [0 0];
mstruct = defaultm(mstruct);
[x,y]=mfwdtran(mstruct,-77.51,166.4)  % produces [0.092723,-1.7802]

ArcMap 10.3.1 code (default parameters, arcPy code)

    arcpy.Project_management("point.txt","PROJECT", "GEOGCS['GCS_Sphere',DATUM['D_Sphere',SPHEROID['Sphere',6371000.0,0.0]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Azimuthal_Equidistant'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")


texts: Why does OWSLib's getfeature(..., filter=...) not filter?

texts:
I am struggling to make a filter effective for a GetFeature request with OWSLib. Here's the code I wrote so far (I'm extracting data from Berlin's open geodata-server):

from owslib.wfs import WebFeatureService
from xml.etree import ElementTree as ET

wfs = WebFeatureService("http://fbinter.stadt-berlin.de/fb/wfs/geometry/senstadt/re_vms_detailnetz_str", "1.1.0")
filter_ = PropertyIsLike(propertyname='BEZIRK', literal='Mitte', wildCard='*')
filterxml = ET.tostring(filter_.toXML()).decode("utf-8")
res = list(ET.fromstring(wfs.getfeature(typename="fis:re_vms_detailnetz_str", propertyname="BEZIRK", filter=filterxml, maxfeatures=10).read()))

This code is copied from the OWSLib sample section on WFS. However, this has not the intended effect, since I get:

print list(list(res[1])[0])[4].tag  # -> '{http://www.berlin.de/broker}BEZIRK'
print list(list(res[1])[0])[4].text # -> 'Spandau' != 'Mitte' -- Oops?


texts: Clustering a number of points with the same coordinate around a central point

texts:
Is there a plugin in QGIS to randomly distribute points with the same coordinates around the coordinate? Similar to the "shift point plugin", but instead of producing a circle of points the points are distributed around the coordinate randomly within a set distance or polygon.

Is there a plugin that does this or is there a combination of tools I could use to do this? I have looked online but can't seem to find anything.


texts: Getting all results from Selection made by Attribute or Location using ArcPy?

texts:
Using arcpy I'm trying to get the results of a selection, or what is currently selected as a list.

Basically I'm attempting to build an addin which on click will either perform an operation where I invoke SelectLayerByAttribute_management (or ByLocation) to select some features, or the user performs a selection ByAttribute (or ByLocation), and while they have their selection hightlighted executes my snippet by pressing a button on the toolbar.  In each  case I want to get a list of results to perform a specific task on only certain selected features.

I can't find out how to get a list of the currently selected features, or find a listener for  selection having been made so I can enable the button, and ensure I can perform my operation on an existing selection.


texts: Get a sql view to persist in a CartoDB visualization

texts:
I'm new to CartoDB and so far impressed. But I'm struggling with how to get a sql query to persist in a visualization. I type it in, like 'SELECT * FROM env_feedlots where active = 'Y' and BIRDS > 0', hit "Apply", and see the effects in the map. But I get a throbber in my browser tab that never seems to die.

And then, when returning to the visualization, it seems the sql query is gone.

Here's a link to my viz: http://cdb.io/1ncUdG7

When you load it up, if you see my sql query, you'll see clusters of dark blue in the central and southeast portions of the state, and very light blues in the southwest. It seems like that works through the "public map" link. But why doesn't it "keep" when I return to the visualization through my dashboard?

Edit:
I found a "fix", if I click "clear this sql view", I can start over, re-type the SQL I want, click "Apply", and it will persist. But that begs the question: why can't I ADD to an existing SQL query and have that addition persist?


texts: Incorrect positioning in Google Earth Engine after exporting from Google Earth Engine

texts:
I have an image loaded in Google Earth Engine and I want to save an enlarged area around it both as a feature collection and another image. I use below code for this purpose and it works fine and objects overlay correctly when I add the image (original image), box_fc, and box2X layers on GEE map. The original image projection is 'Albers_Conical_Equal_Area'.

var proj = image.projection()
var box = image.geometry()
var features = [ee.Feature(box.buffer(2500).bounds(0.001, proj=proj), {'Value': 0})];
var box_fc = ee.FeatureCollection(features)
var box2X = box_fc.reduceToImage({
  properties: ['Value'],
  reducer:ee.Reducer.first()
  }).reproject(proj)
Export.table.toDrive({
  collection: box_fc,
  description: 'box2X',
  fileFormat: 'SHP'
})
Export.image.toDrive({
  image: box2X,
  description: 'box2X',
  region: box_fc,
  scale: 30
})

However, when I load the exported shape and image I see they are far apart. The shapefile (box_fc) is aligned correctly with the original image, but image (box2X) not. When I look at the image properties, the projection and extent are also correct, but it fails to display it correctly.

What is wrong here?

On the shapefile properties in ArcGIS Desktop, there are two question marks (??) after extent values.

What is it for?

I found that actually the 'box2X' image has not been exported properly and it looks like my understanding of reprojection is not correct. When I load the exported 'box2X' again into Google Earth engine and display it on the map, it is displaced similar to what I see in ArcGIS Desktop. I don't understand why before export it is displayed correctly in GEE, but not after export.


texts: Creating Thiessen (Voronoi) polygons using lines (rather than points) as the input features?

texts:
I have a set of line features inside a particular polygonal boundary. For each line, I'd like to generate a polygon inside which every possible point is closer to the given line than to any other line in the layer. I've done this in the past for point input features using Delaunay triangulation, but if there's a similar process for doing it with line features I haven't been able to find it.

ETA: Geogeek's solution had occurred to me, but in straighter sections where the input lines  have fewer vertices, the resulting polygons get way too close (even overlapping) a line that they shouldn't. Here, the red lines are my inputs, you can see the vertices and the Thiessen polygons generated from them.

Perhaps a quick and (very) dirty solution might be to convert each line to a plentiful set of evenly-spaced points (rather than the line's vertices only), generate Thiessen polygons from those, then dissolve them based on the originating line ID.


texts: How to determine aerosol model value for i.atcorr in GRASS? (Using Aeronet data)

texts: I am trying to conduct atmospheric correction of landsat 8 image using i.atcorr of Grass GIS. In thıs regard, I need to know Aerosol Optical Depth (AOD) as stated in Grass wiki page. I've found some valuable leads ın this site and explored Aeronet (http://aeronet.gsfc.nasa.gov/new_web/aerosols.html). I prefer it (relative to MODIS MOD08) since it offers a lot more detail and timeseries of data (I think I just need the data for the time when image was taken). However, the data file for my specific site contains no data at AOT_550 (I presume that is the exact wavelength sought). there is AOT_555 and AOT_551 but without any data. The closest data I found is AOT_500 and AOT_675. Another thing I noticed, the data value gradually decreases from AOT_340 to AOT_1640. My question, Is it possible to use a linearly interpolated value for AOT_550 by using AOT_500 and AOT_675 data?


texts:
I am getting GeoJSON data from postgresql, but it seems to be invalid. How can I validate it?

After getting GeoJSON format data I want to put on map for projection. I am using LeafLet, java EE, postgresql.

<%@page import="java.io.Console"%>
<%@page import="org.json.simple.JSONArray"%>
<%@page import="org.json.simple.JSONObject"%>
<%@page import="java.util.ArrayList"%>
<%@ page language="java" contentType="text/html; charset=ISO-8859-1"
    pageEncoding="ISO-8859-1"%>
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">

<meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1">
<title>Insert title here</title>

<%
JSONObject GJson = (JSONObject) request.getAttribute("geoJson");
//String GJson = (String)request.getAttribute("geoJson");
System.out.println("data at jsp : "+GJson);
%>

GeoJSON retriving code(edited):

<%
//JSONObject GJson = (JSONObject) request.getAttribute("geoJson");
ArrayList GJson = (ArrayList)request.getAttribute("geoJson");
//String GJson = (String) request.getAttribute("geoJson");
//String GJson = (String)request.getAttribute("geoJson");
System.out.println("data at jsp : "+GJson);
%>


texts:
I just started using google map API but found some restrictions on daily request limit. After some research I found I can use open street map, but has  its own request restrictions.

To resolve I installed nominatim on my server and started using that until I found out in non metro area nominatim has very few geo points but available in google map.

I am thinking to use below approach:

Setup Nominatim on my server.
Insert OSM data to database.
Sync database with nightly OSM data.
Take the geo point from google and insert into my database.

I have below queries regarding this approach:

Is there any nominatim API which I can use to insert those new geo points from google in my nominatim db.
How will I resolve conflicts between OSM and google as I will be constantly syncing with OSM data.


texts:
I want to cache my raster data(large geotiff files) on local database. I thought to use spatialite. I found librasterlite which provide functionality to do that.

Can anybody point to some tutorial on it.
Or is there any other alternative.

Programming environment - C++

I want to store geotiff data to spatialite database. I am using gdal to read metadata and reading tiles converting them to hexwkb format and storing it in raster table (id (integer) and raster (blob)  as columns).Table is created but when i am opening my database and select raster data column it shows nothing. I mean i can not explore data and see the tile image. Can any body tell me what should be the format in which we should store the tiled data. Are there any function which can automate these things. Like I tried to integrate rasterlite_load in my code but there were memory leakage in the code. Is there any function in spatialite or rasterlite which can take char buffer of tiled image and store it into blob data type in table (and corresponding function for reading).


texts: How to amend Python code to work in QGIS 3?

texts:
In another thread, I was given this code to import large numbers of dxfs into QGIS. This code works in QGIS2.18, but they just released QGIS 3.0 and have changed the codes to use a different version of Python. This code doesn't work in 3.0. Running it generates a series of requests to set the CRS, but nothing then imports.

How can I edit this code (or other codes) to be compatible with QGIS 3.0? What was changed?
I tried an online converter but it didn't alter anything.

import glob, os
path = "E:/aawork/aasites/aa-all towns and villages/Hemelhempstead/osmap 2012/"
for layer in glob.glob(path + "*.dxf"):
    vlayer = QgsVectorLayer(layer, 'name', 'ogr')
    subLayers = vlayer.dataProvider().subLayers()
    for subLayer in subLayers:
        geom_type = subLayer.split(':')[-1]
        if geom_type == 'LineString':
             uri = "%s|layername=entities|geometrytype=%s" % (layer, geom_type,)
             dfx_file_name = os.path.splitext(os.path.basename(layer))[0]
             layer_name = "%s - %s" % (dfx_file_name,geom_type,)
             sub_vlayer = QgsVectorLayer(uri, layer_name, 'ogr')
             QgsMapLayerRegistry.instance().addMapLayer(sub_vlayer)


texts: Why does my WMS layer overlap the Google Maps v3 logo?

texts:
I'm using Google Maps v3 as a base layer and on top of that I create an overlay pointing to my WMS (streets and avenues).

A few months ago all works fine, my WMS layer was not overlapping the Google Maps footer (Google logo and copyright), but today it is.

Here is how it looks:

Here is the code:

var
  map,
  stylez = [
    {
      featureType: "administrative",
      elementType: "geometry",
      stylers: [
        { visibility: "off" }
      ]
    },{
      featureType: "landscape",
      elementType: "geometry",
      stylers: [
        { visibility: "off" }
      ]
    },{
      featureType: "road",
      elementType: "geometry",
      stylers: [
        { visibility: "off" }
      ]
    },{
      featureType: "transit",
      elementType: "geometry",
      stylers: [
        { visibility: "off" }
      ]
    },{
      elementType: "labels",
      stylers: [
        { visibility: "off" }
      ]
    }
  ],
  gmap = new OpenLayers.Layer.Google("Google Styled", {
    type: 'styled'
  }),
  styledMapOptions = {
    name: "Styled Map"
  },
  styledMapType = new google.maps.StyledMapType(stylez, styledMapOptions),
  projection = new OpenLayers.Projection("EPSG:900913"),
  displayProjection = new OpenLayers.Projection("EPSG:4326"),
  bounds = new OpenLayers.Bounds(-57.81887, -25.53276, -57.24308, -25.17649).transform(
    displayProjection,
    projection
  );

// avoid pink error tiles
OpenLayers.IMAGE_RELOAD_ATTEMPTS = 3;
OpenLayers.Util.onImageLoadErrorColor = "transparent";
OpenLayers.ImgPath = "http://js.mapbox.com/theme/dark/";

function init() {

  var
    options = {
      controls: [
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.PanZoomBar(),
        new OpenLayers.Control.LayerSwitcher(),
        new OpenLayers.Control.MouseDefaults(),
        new OpenLayers.Control.KeyboardDefaults()
      ],
      projection: projection,
      displayProjection: displayProjection,
      units: "meters",
      numZoomLevels: 12 /* 18 */,
      minZoomLevel: 6,
      maxResolution: 156543.033928041,
      maxExtent: new OpenLayers.Bounds(-2.003750834E7,-2.003750834E7,2.0037508345578495E7,2.0037508345578495E7)
    };

  map = new OpenLayers.Map('map', options);

  // base layer
  var tiled = new OpenLayers.Layer.WMS(
    "Roadmap",
    "http://foo.com:80/geoserver/gwc/service/wms",
    {
      LAYERS: 'Foo',
      STYLES: '',
      format: 'image/png',
      tiled: true,
      transparent: "true",
      tilesOrigin: map.maxExtent.left + ',' + map.maxExtent.bottom
    },
    {
      isBaseLayer: false,
      opacity: 1.0,
      visibility: true
    }
  );

  tiled.events.register("loadstart", tiled, function() {
      console.debug("Foo, Load Start");
  });

  tiled.events.register("loadend", tiled, function() {
      console.debug("Foo, Load End");
  });

  map.addLayers([gmap, tiled]);

  map.zoomToExtent(bounds);

  gmap.mapObject.mapTypes.set('styled', styledMapType);
  gmap.mapObject.setMapTypeId('styled');

}


texts: Great-circle lines in Equirectangular projection

texts:
Just to check that I'm on the right track:

Are all great-circles on the sphere and in equirectangular projection (i.e. latitude, longitude pairs) either:

meridians (i.e. going pole to pole)
of the form tan latitude = sin360(longitude + rotation) * amplitude + offset

(with additional restrictions on the offset/amplitude combinations -- obviously, all great circle paths with an amplitude of 0 also have offset 0 - the equator).

Or are there great-circle paths that do not fit into this scheme (again, only in a longitude-latitude coordinate system, not on other map projections).

Note: I added the tan above after posting the question, in a reaction to whubers excellent reply. It turns out that offset then always is 0.


texts: Spatial analysis of classified grid using QGIS

texts:
I'm calculation the positional accuracy of linear features using the buffer method and dividing the region using a rectangular grid.
Here follows the vectos to be compared.

The grid cells were classified according to difference of vector lengths inside the cell, so the red cells have much more of one vector than another, blue cells are the opposite and yellow cells have similar lengths of both vectors.

I'm testing different cell sizes for the grid.
The smaller the cell size is, the better will be the "resolution" for the spatial analysis, so will be possible making inferences of more detailed objects, but when I use a too small cell size, the matrix looks like the vectors themselves, making useless using a matrix instead of vectors.

Here is the grid with cell size 0,001 degree.

Here is the grid with cell size 0,0025 degree.

Here is the grid with cell size 0,005 degree.

Here is the grid with cell size 0,0075 degree.

So there seems to have a optimum size of cell, which is the smallest that don't "look like" the vectors that I want to determine.
In this case is 0,005 degree (which was the same size when I tested with some other datasets).

My question is: There is a way to measure that the grid "look like" the vectors? Because currently I'm just using my personal perception instead of science.

Just need a scientific way to show that the 0,001 and 0,0025 degree cell grids looks like the vectors and the 0,005 and 0,0075 degree cell grids doesn't. I can't write in an article paper that the grid "looks" like the vectors because it's too subjective, I need numbers show that.


texts: Should 90° S, 180° W and 90.0°S, 180° E produce different C-Squares

texts:
I created a simple C-Squares encoder in Python as a training exercise to help me understand how the encoding system works.

From my understanding limited understanding of geosystems and reading the spec;
    90° S, 180° W and 90.0°S, 180° E are both the same spot on the earth and yet the C-Squares are '5817:499' and '3817:499' respectively.

Is this correct?

https://github.com/zeonglow/csquares


texts:
When i create widgets for the attribute form in the properties of a layer:

Is this stored in the geopackage layer or in the QGIS project file?


texts: No cadastral maps in the UK. Solutions/alternatives for a methodology?

texts:
I need help reworking a methodology in the absence of cadastral data in the UK.

The aim is to plot a range of businesses in a UK city on separate layers (according to business type) as polygon shapefiles for spatial analysis down the track.

The process would go from a database created or search results saved as KML in Google Earth and then uploaded to either QGIS or ArcGIS. From there I would need to create a polygon around the complete building footprint of that business premises.

I need a relatively fast system as there may end up being 1000's of business to plot.

If I had cadastral maps this would be much faster in terms of tracing etc. Alas, the UK doesn't have any official cadastral data.

Any ideas on how to refine this methodology?


texts: Move cached layer to another geoserver instance

texts: I have cached a layer on a local geoserver using the internal GWC. I want to publish the created cache on another geoserver for production. I have moved the tiles but geoserver does not pick it up the added cache. I assume I need to move or create/modify some configuration for geoserver to be aware if it but I am not sure what I need to do?


texts:
I use this code in order to remove layers that are out of the data frame of the map:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.overwriteOutput = True
env.workspace = r"G:\PROJECTS\kohav_yair_west\gis\aa"
for mxdname in arcpy.ListFiles('*.mxd'):
        print '---' ,mxdname, '---'
        mxd = arcpy.mapping.MapDocument(os.path.join(env.workspace, mxdname))
        df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
        print 'REMOVED LYR'
        for lyr in arcpy.mapping.ListLayers(mxd, "" ,df):
                if lyr.isGroupLayer == True: continue
                if lyr.isRasterLayer == True: continue
                arcpy.MakeFeatureLayer_management(lyr.dataSource, "lyr")
                frameLyr = r"G:\PROJECTS\kohav_yair_west\gis\layers\frame.shp"
            arcpy.SelectLayerByLocation_management("lyr", "INTERSECT", frameLyr)
                if int(arcpy.GetCount_management("lyr").getOutput(0)) == 0:
                        arcpy.mapping.RemoveLayer(df, lyr)
                        print u'Removed '+unicode(lyr) + ' lyr'
        print '---REMAIN LYR---'
        for lyr in arcpy.mapping.ListLayers(mxd, "" , df):
                print lyr.name
        mxd.save()
del mxd

but i get this error and i don't understand what is the problem?

--- ttl.mxd ---
REMOVED LYR

Traceback (most recent call last):
  File "G:\desktop\y\yaron\Private\shonot\software\---gis---\tools\YARON_SCRIPTS\SelectLayerByLocation3.py", line 21, in <module>
arcpy.MakeFeatureLayer_management(lyr.dataSource, "lyr")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 6520, in MakeFeatureLayer
raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset F:\GIS\Updating_2014\ttl\arcview\ttl_all_blue_lines_07.05.shp does not exist or is not supported
Failed to execute (MakeFeatureLayer).

This route : F:\GIS\Updating_2014\ttl\arcview\ttl_all_blue_lines_07.05.shp exist but the name of the layer is :"ttl_all_blue_lines_07.05.14.shp" and in the error it mention as ttl_all_blue_lines_07.05.shp


texts: Extracting Summary Statistics value with cursor and assigning it to a variable

texts:
I am writing a script (to be used in ArcGIS 10.7.1) that uses a Zonal Statistics tool to create an average of a raster that is beneath a grid. I then use the Summary Statistics tool to take the summation of all the grid values and display them in a table. The resulting table has only one row with three columns, as shown below:

I'm trying to use a cursor to assign the value of the "SUM_Value" to a variable, because that value will be used later to compute some fertilizer amounts.

How can I do this?

The methods that I've tried aren't working out. I get this error when I do so:

Traceback (most recent call last):
  File "<pyshell#132>", line 1, in <module>
    totalNitrogen = row.getValue(field)
NameError: name 'row' is not defined

Here is what I've tried:

#Define workspace
env.workspace = "C:\\someWorkspace"

#Defining variable 'Grid'
Grid = #Previously created grid

#Zonal statistics tool
recommendationMap = ZonalStatistics(Grid, "grid_field", outputRaster, "MEAN")

#Compute raster to Int
recommendationMapInt = Int(recommendationMap)

#Save raster
recommendationMapInt.save(env.workspace + "\\nitrogenRecommendationMap")

#Generate statistics with Summary Statistics tool
nitrogenStats = arcpy.Statistics_analysis(recommendationMapInt,"outputRaster",[["Field_to_summarize", "SUM"]])

#Resulting output field defined for cursor
nitrogenField = "SUM_Value"

#Cursor defined
nitrogenCursor = arcpy.SearchCursor(nitrogenStats)

#Error occurs here
totalNitrogen = row.getValue(field)

As previously stated, I need that totalNitrogen variable to equal the "SUM_Value" field from the Statistics_analysis table.


texts: Convert a raster with one band into a raster with RVB bands on QGIS?

texts: So yesterday i downloaded a raster tile on a WMS server and exported it with Grass. I get a raster with 3 bands, i treat with the Semi automatic classification Plugin in QGIS and it works perfectly. But today i realized that said raster didn't cover the whole area i'm working on right now, i download the raster again, from the SAME source with the SAME settings (i checked twice) i export it using GRASS again and when i try to resample said raster with SCP it doesn't recognize it. I don't uderstand why imediatly up until i check the bands and realize that now, the raster is on only one classified band and not on RVB. Since the raster took hours to download i would like to know if i could convert this panchromatic tile into an RVB one and, if their is no simple way, how to force the r.in.wms function to take RVB rasters when it is avaible.


texts: Georeferencing JPEG from TIFF using ArcMap?

texts:
I recently georefernced some TIFFs. However, unfortunately I have to do the same stuff with equivalent JPGs (same image, same resolution, just JPGs).

Is there an easy way to georeference the JPGs using the TIFFs using ArcMap 10.3?


texts: Subtracting the current NDVI image from the previous NDVI temporal cube image

texts:
I am trying to perform change detection by subtracting indexes.

I am still new to GEE.

I am taking advantage of the structure of an existing script to try to achieve my goal.

I am generating an image from NDVI-quality temporal cube (qualityMosaic).

I cannot subtract each 'current' NDVI-image from the previous NDVI-image temporal cube image.

How can I concatenate this subtraction?

I forgot to say that I need the result to amaze the date metadata.

Below is the script.

// Define time interval to composite scenes
var date_start="2018-01-01";
var date_end="2018-12-30";

// define mask functions

function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}
// empiric mask
var maskNotValid = function(image) {
  var cloudMask = image.select(Band_BLUE).lt(1500).and(image.select(Band_LSWIR).gt(600))
  return image.updateMask(cloudMask);
}

var L8coll = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filterBounds(Region)
.filterDate(date_start,date_end)
.map(maskL8sr)
.map(function(image){
  return image.rename(['blueblue', 'B1', 'green', 'B3', 'B4', 'B5', 'B7', 'brigtemp', 'B8', 'B9', 'pixel_qa', 'B11']);//padroniza a B7
})
.select(['B1', 'B3', 'B4', 'B5', 'B7', 'pixel_qa']);

// merge collections
var image_collection_type = L8coll;
print (image_collection_type)

// define band names
if (image_collection_type){
  var Band_BLUE='B1'
  var Band_RED='B3'
  var Band_NIR='B4'
  var Band_SSWIR='B5'
  var Band_LSWIR='B7'
  var resolution=30
}

// Calc spectral index
var add_NDVI = function(image) {
  return image.addBands(image.expression('float(b("'+ Band_NIR + '") - b("' + Band_RED + '")) / (b("' +  Band_NIR + '") + b("' + Band_RED+ '"))').rename('NDVI'))
};

//Add date variable
var add_date = function(image) {
    var timestamp = image.metadata('system:time_start');
    return image.addBands(timestamp);
};

var get_date=function(d){
    return d.set("date",ee.Date(d.get('median')).format('YYYY-MM-dd'))
};

// Change Detection
var addChangeQualityBand = function(image) {
    return image.addBands(image.normalizedDifference([Band_LSWIR, Band_NIR]).rename('Change'))
};

var addNotChangeQualityBand = function(image) {
    return image.addBands(image.normalizedDifference([Band_NIR, Band_RED]).rename('NotChamge'))
};

// Define Maximum NBR
function get_post_composite (date_start, date_end){
  var image_collection = ee.ImageCollection(image_collection_type)
    .filterDate(date_start, date_end)
    .map(maskNotValid)
    .map(addChangeQualityBand)
    .map(add_NDVI)
    .map(add_date)
    .filterBounds(Region);
  var composite = image_collection.qualityMosaic('Change');

  return composite
}

// Define Maximum NDVI
function get_pre_composite (date_start, date_end,image_collection_type){
  var image_collection = ee.ImageCollection(image_collection_type)
    .filterDate(date_start, date_end)
    .map(maskNotValid)
    .map(addNotChangeQualityBand)
    .map(add_NDVI)
    .map(add_date)
    .filterBounds(Region);
  var composite = image_collection.qualityMosaic('NotChamge');
  return composite
}

var post_composite=get_post_composite (date_start, date_end,image_collection_type)
var pre_composite=get_pre_composite (date_start, date_end,image_collection_type)
var subtract_post_pre_composite=post_composite.subtract(pre_composite)

//clip
var post_composite = post_composite.clipToCollection(Region)
var pre_composite = pre_composite.clipToCollection(Region)

//plot
Map.addLayer( pre_composite , {bands: [Band_LSWIR, Band_NIR, Band_RED], min: [0, 0, 0], max: [4000, 4000, 4000]}, 'NDVI maximum composite');
Map.addLayer( post_composite , {bands: [Band_LSWIR, Band_NIR, Band_RED], min: [0, 0, 0], max: [4000, 4000, 4000]}, 'NBR maximum composite');


texts: I am stuck with a task to compare two line features from 2 different shape files. They don't have any common field so i can't join them. I have source for street names and existing street names in database. My task is to compare which of these names match and find the mismatches.


texts: How to set NULL as default value for a date field in QGIS 2.4?

texts: In QGIS 2.4, i am unable to set the default value as NULL if i have not updated any value on date fields. Its automatically taking the current system date. I am expecting to set the default values(NULL) for date fields in properties of layer, since that default values should be displayed while capturing the new feature itself. It may be edited if required. Please refer the attached screen shots.


texts: Using ArcPy and ArcGIS API for Python together?

texts: I am having trouble with how I am supposed to use both ArcPy and ArcGIS API for Python. ArcPy is only in Python 2 and ArcGIS API for Python is in ArcGIS Pro and is Python 3. I don't know how to set my paths in VS Code or Atom to be able to use both.


texts: Move GeoServer folder inside Jetty webapps to Tomcat webapps folder

texts: Can I move the GeoServer folder app directory to the Tomcat webapps directory to migrate it to Tomcat? I do not want to update GeoServer version or do any configuration change, only migrate it to Tomcat.


texts:
I'm trying to test an mxd with mxdperfstat.

I want to test several scales, but as soon as I provide a scale over 500, the program simply quits.
No Error warning. It just stops working.
If I provide several scales at once (250;500;800;1000) it tests the 250 and 500 and then stops.

I use ArcGIS 10.2.2 and the according mxdperfstat10.2.exe.
Thats how I call the programm:
mxdperfstat10.2.exe -mxd F:\Path\test.mxd -scale 250;500;800;1000 -xy 4348929;5498749 -width 1000 -height 1000

The end of the console output looks like this:

21.11.2014 13:25:18 RunTest:LoadXML
21.11.2014 13:25:18 RunTest: ESRI.ArcGIS.ADF.COMSupport.AOUninitialize.Shutdown()
21.11.2014 13:25:18 RunTest: ESRI.ArcGIS.ADF.COMSupport.AOUninitialize.Shutdown() Completed
21.11.2014 13:25:18 RunTest:Application.Exit()
21.11.2014 13:25:18 After RunTest:Application.Exit()

No error is documented.

In the resulting .xml report there are only the test values for scale 250 and 500.
The behavior is the same for different MXDs.

Any Idea what the Problem could be here?


texts: Windows Embedded Handheld 6.5 Professional - How to Uninstall VBScript Runtime Library and Install again?

texts:
I've been having trouble getting the Trimble Positions for ArcPad Extension to correctly install on my Trimble Geo 7X handheld. I saw that another user asked about something similar and the solution provided was to re-install VBScript Runtime Library. I've gotten that far, but can't find a way to re-install it.

What is says:

Script Error 800A01A8

Error Source Unavailable

Error Description Unavailable

[Line: 61, Column: 4]

Source Text Unavailable

Also, upon exiting ArcPad, it throws another error popup that references Line: 41, Column: 5. Has anyone experienced anything similar?


texts: Apply symbology to layer using arcpy - emulating the "import symbology function"

texts:
I am using ArcMap 10.4.

I am creating a Python Toolbox which runs through a series of calculations and geoprocessing tools. The Toolbox is an update of an old VBA script I ran in ArcMap 9.3, where it ran perfectly. I am updating it to run in ArcMap and Pro.

One of the functions used is intersecting a user defined polygon with a larger base layer. This base layer contains a number of fields, including a code which I use to apply a Symbology to it.

I am having an issue assigning the symbology of the larger Polygon Layer to the smaller intersected Polygon using the arcpy.ApplySymbologyFromLayer_management statement.

See an example of the base layer below:

Running the Toolbox produces this result:

The blue area should be the speckled the white and brown from the first picture.
I thought maybe there was an issue with the attribute code I was using to assign the symbology, but when I used the "Import Symbology" under the Symbology Tab, the correct symbology was displayed.

The attribute for the symbology is a 5 digit number, it is static, so will never change.

The symbology itself is the "Unique Value" option.
The code I am using is below. I have already imported arcpy and completed the intersection before this code runs.

    #File to symbolise
    Sym_In = FileN
    #Location of the lyr file to use
    Sym_Lay = r"C:\Data\Layer_v2.lyr"
    #Even if I put the field containing the code in here...nothing changes
    Sym_Field = ""
    arcpy.ApplySymbologyFromLayer_management (Sym_In , Sym_Lay)

Is there a way to get the result of the applysymbologyfromlayer to match the "Import Symbology" button in the layer properties?


texts: formula to calculate geo-referencing paramerters

texts:
I am trying to overlap the vector and raster data. I have two control points on both raster and vector data. Now I want to generate world file parameters for overlapping two data. I am using DotSpatial (github DotSpatial) as an base library for rendering the GIS data.

My control point on raster image are: (x1,y1) and (x2,y2) and control point on vector data are: (u1,v1) and (u2,v2).

How can I generate the required parameters programmatically?

Parameter details are here

I need something like below


texts: How to setup Visual Studio Code (Windows) for PyQGIS 2 or 3?

texts:
I'm trying to use Visual Studio Code (VSC) as IDE for my QGIS plugins with the Python extension, but something doesn't work as expected. I've been reading about the PyCharm config from Gary S. but I couln't apply it to VSC.

Should I use an .env file? What should I set in the settings.json? How to call batch files when VSC starts?


texts: Why do St. Croix and Angaur of Virgin Islands have same FIPS code?

texts:
Could any one help me answer this.
I am a non US citizen and working with d3js(geojson). I came across two locations in Virgin Islands which share the same FIPS code - St. Croix and Angaur have FIPS code 78010. I am looking for some sort of explanation around this.

I couldn't find any answer on the web.


texts: Labels not displaying after calling arcpy.da.SearchCursor

texts:
I'm attempting to create labels for a feature class that labels each feature "[ParentFeatureName] - [ThisFeatureID]". I found a code snippet on Esri support that should fit the bill (with the variable names and return changed as needed):

def FindLabel ([keyField], [FirstLabel]):
    import arcpy

    key1 = [keyField] # Key field in feature class
    key2 = "ID"   # Key field in related table
    L = [FirstLabel] # Label field in feature class
    L2 = "Label2"   # Label field in related table
    myDataTable = r"ExampleTable"   # Path to related table

    cur = arcpy.da.SearchCursor(myDataTable, [key2, L2])
    for row in cur:
        if str(key1) == str(row[0]):
            L = L + " " + str(row[1])
    return L

But no labels appeared! So I did some experimenting:

I did a sanity-check in the Python pane to make sure that the comparison is valid (it is)
I started commenting out code to see where it bombed.

I had expected that commenting out the for loop and the assignment of L would label my features with L, but to my surprise, there were still no labels. Only when I also commented out cur = arcpy.da.SearchCursor(myDataTable, [key2, L2]) did the labels appear!

The same behavior also happens using the legacy arcpy.SearchCursor.

Why is that line of code preventing the labels from appearing and how do I work around it?


texts:
I have dump my data in PostGIS from shapefile successfully. I made a postgres connection through PHP and that also returning the data of other field. But the main problem is GEOM field visualization. I want to visualize GEOM field (Polygon and Polyline). Simply I want to visualize PostGIS data to web. Any kind of help or tutorial will be appreciated highly.

I am clarifying little more. I have created my own tiles with tilemill. I am now using Leaflet to visualize the tiles. I have completed this portion successfully. Now I have same database in PostGIS.

I want some interactivity with that data. Such as Search return polygon/polyline/point on the map from PostGIS database. What is the easiest way to do this.
In future I want to integrate PGRouting with this map also.


texts: Automated work flow web form - QGIS Spatial Query - PDF report output

texts:
With some limited QGIS experience is it possible to create an automated workflow that takes a location from a form on our website into QGIS, does a couple of spatial queries then kicks out an image and CSV table into a pdf report that can then be emailed back to a user of the website?

I'm working on macs so would rather not get into ESRI products.

Some initial research leads me to think that a combination of Applescript, QGIS and Filemaker may be able to achieve this however I would also welcome any suggestions for software solutions like FME.

Apologies if this question is too general for the forum.


texts: Editing values in attribute table column

texts:
I have an attribute table with 5,000 records in QGIS. From this I have selected a subsection of 200 records whose attributes I need to modify. In the column called 'NAME' I would like to replace the current values with the same value for all 200 selected records.

Right now, the values in column NAME are different for each row (eg. Tom, Claire, Michael etc.) and I would like to know if there is a way to replace the 200 random names with the same name. I don't want to create an additional column, only to overwrite the current value stored in the 200 selected records.

So far I've tried various replace and regexp_replace expressions but I can't seem to get the correct syntax. Is there syntax that is able to select all the random names in the subset and replace it with a name I have chosen. Something like select, from selected NAME records and replace with XYZ?


texts:
It seems that the OpenLayers.mobile.js library (v 2.12) doesn't have a popup class. Does anyone know how to do a popup with the mobile version?

I'm using some code that i ported over form normal OL web map. When i call new OpenLayers.Popup.Anchored i get an undefined error on OpenLayers.Popup class.

It's not in the OpenLayers.mobile.js. I'm lost


texts: How to put North Arrow and Scale on map in COMPOSER in QGIS 2.2.0

texts:
I am new to QGIS and the questions asked and answered here about the North Arrow do not seem to help me solve my problem below.

In the main mapping pane of QGIS 2.2.0 Valmiera, I can add the scale and north arrow, but as soon as I go to the Composer to print a map, those two things disappear from the map that I open in the composer.

The Composer has something which allows me to add a scale thing, but it shows the actual scale on the paper that I'll be printing on and not the actual scale of the surface of the earth of the displayed map. I even tried refreshing the map but it did not work.

The bottom line is this:  How do I get a scale thing (sorry!) (expanding and contracting in values as I zoom in and out) in Composer, and how do I get a North Arrow in there as well.


texts: Leaflet Routing Machine show option doesn't work

texts:
I am adding route to the map and I want to hide instruction (turn left, turn right etc.) but I want to show them later. As I can see there is parameter show for that. If I set show: false I see white squares instead of istructions which is ok, the problem is that I can't show instructions then by clicking on that squares, no interaction with them at all.

L.Routing.control({
  waypoints: [
    L.latLng(54.736985, 55.980878),
    L.latLng(54.732798, 55.969934),
  ],
  show: false,
}).addTo(map)

There is nothing wrong here I guess. What can it be then?


texts: Map atlas from multiple layers but same extent

texts: I need to create and export several maps, one from each raster layer in the TOC (~50, so preferably not manually) and I want the extent to remain the same. I have explored map series function with ArcGIS Pro but cannot figure out how to do it, as it only seem to work with different zoom/extents showing the same data layers. Any idea on how to accomplish this? Happy to try if some coding is required but please provide some hints as I am a beginner with ArcPy.


texts: Are there any well symbols for ArcGIS Silverlight/WPF?

texts: I'm looking for some standard map symbols for a silverlight app.  Does anyone know where I could find these?


texts: Splitting all lines in layer at vertices using QGIS

texts:
I want to automatically split all lines in a QGIS layer at their vertices (like this).
I've done a Google search and found this article, which is quite close but the workflow is complicated. Because I'll have to train my co-workers (who are unfamiliar with GIS) to use this tool, I want something as simple as a plugin or something that doesn't require GRASS.
Since it's been a long time from 2010 - the time of the previous answer, I wonder if there is any better way to achieve the same result?


texts: How to debug Live GPS and base layer alignment issues in QGIS?

texts:
I am trying to use a bluetooth GPS with a couple of base layers in QGIS. The GPS/bluetooth connection seems to be working great...I can connect and disconnect without issue. I've set my project CRS to WGS84, and I've enabled OTF CRS transformation. I've confirmed the CRS for my base layers is set properly. However, I can't zoom to my base layers to know if my GPS point is lining up where I think it should be with the base layers. Whenever I click 'Zoom to Layer Extent', nothing happens. Any thoughts on what I might be doing wrong?

I am using QGIS ver. 2.0.1 on a Windows 8 32-bit tablet.


texts:
I have GPS tracking solution and most my customers in US. My reverse-geocoding solution pretty basic:

Table with US ZIP database
Procedure in SQL Server where I search for city/state closest to reported lat/lon

Now I need to add Canada to my service coverage. I will buy Canadian postal codes and plan to store it in separate table.

So, straight-forward solution will be to:

Check if there is close US city/state and report
If not found - check in Canada table

This solution will put extra load on server for those canadian users since US table will be queried first.

For various reasons I don't want to merge tables.

However, if you look at map - you can tell that Canada is above USA :)

So, I was wondering if I can use some hard-coded numbers to know immedeately where this point is.

Right? Wrong? Other ideas on how to tackle this?


texts: Why does the Translation from 3 Band TIFF to ECW in FME 2017 fail?

texts:
I have a number of 16bit RGB raster images which need to be converted to ECW format. The images are displayed perfectly in ArcGIS and have the following attributes.

rows and columns: 1506, 2051
No Bands: 3
Format: TIFF
Pixel depth: 16 Bit
Compression LZW

When I load the Tiff in FME 2017 and inspect it, it is rendered black with the following information:

rows and columns: 959, 250
BAND 0 (Gray 16)
BAND 1 (UINT 16)
BAND 2 (UINT16)

My translation to ecw fails due to "Unsupported band interpretation 'GRAY16'".
How can I resolve this?


texts: No sync geopackage raster layer to QFIELD

texts:
I have a problem with exporting raster layers stored in a geopackage, to a QFIELD project. I work with version 3.4.12 of QGIS and I have version 3.5.1 of QFIELD Sync installed.

The process I follow is as follows:

I convert a TIF image into a layer of a geopackage

gdal_translate --config OGR_SQLITE_SYNCHRONOUS OFF -co  APPEND_SUBDATASET=YES -co RASTER_TABLE:foto -co TILE_FORMAT=PNG_JPEG -tr 0.5 0.5 -a_srs EPSG:25829 -of GPKG foto.tif geopackage01.gpkg

I build the pyramids of the image

gdaladdo --config OGR_SQLITE_SYNCHRONOUS OFF -r AVERAGE geopackage01.gpkg:foto  2 4 8 16 32 64 128 256

I save the QGIS project, and activate the QFIELD Sync plug-in in the project synchronization configuration options. In the layers list, the one within the geopackage, has the option of no action and does not allow it to be included in the project that will be used by QFIELD.

I think the statements used to create the geopackage are correct, as it does not give any error in QGIS. If I load the original layer in TIF, I can synchronize it without problem. It does not happen to me with a layer only, but it is with any raster that I load from geopackage.

Can anyone have any idea what can happen?


texts: Understanding ArcGIS Flow Direction pixel values?

texts:
The range of values in my flow direction raster  is 1-129.

D8 method was used.

I would expect 128 to be the largest value in accordance with the documentation

what's more I show values OTHER THAN 1, 2, 4, 8, 16, 32, 64, 128

Why are there these other values?

One answer to this question begins to suggest this is normal but does not really explain why/how other values could be produced.
Flow Direction Raster Value Range


texts: Is ArcFM or ArcMap leaking memory when switching between page templates?

texts:
Can anyone else reproduce the memory leak described below? Any insight into what may be causing the leak is appreciated.

I have custom code that swaps between landscape and portrait templates and exports a large number of maps to PDF. Unfortunately, there is a memory leak occurring when I swap from template to template that eventually causes the ArcMap to crash.

I can reproduce the memory leak using the Minerville data and workflow described below.

Configuration:

ArcFM Solution Desktop 10.1

ArcGIS 10.1 for Desktop 10.1

Windows 7 64-bit or Windows XP 32-bit or Windows 2008 64-bit

Workflow:

1) Open Minerville database in ArcMap

2) Open “Stored Items” dropdown

3) Select “Page Templates” tab

4) Double-click the “Adjacent Grids” template to load it

5) Repeat steps 2-3 and double-click “Designer Detail Map Sheet”

6) Repeat steps 2 – 5 and watch memory usage climb in task manager

Using UMDH, the largest ArcFM-related offenders seem to be these:

Stack #1:

ntdll!RtlAllocateHeap+00000274
KERNELBASE!LocalAlloc+0000005F
GDI32!pmfAllocMF+0000018C
GDI32!SetEnhMetaFileBits+00000033
OLEAUT32!_PictLoadUnknownMetaFile+00000088
OLEAUT32!_PictLoadNewImage+000001A0
OLEAUT32!_PictLoadPicture+00000025
OLEAUT32!CPicture::Load+00000027
SystemCore!ObjectStream::LoadObject+000003E9
AfCore!LoadObject+0000005E
RasterGraphicElements!BasePictureElement::Load+000000ED
SystemCore!ObjectStream::LoadObject+000003E9
AfCore!LoadObject+0000005E
AfCore!ExtraSafeReadObject+00000072
PageLayout!PageLayout::Load+00000174
SystemCore!ObjectStream::LoadObject+000003E9
MMStoredStuff!DllGetClassesForCategoryReg+00011B32
MMStoredStuff!???+00000000 : 2DB04186
MMStoredStuff!DllGetClassesForCategoryReg+0000D9E5
MMStoredStuff!DllGetClassesForCategoryReg+0000C297
mmPageTemplatesUI!DllCanUnloadNow+00020402
mmPageTemplatesUI!DllCanUnloadNow+00038EAF
mmFrameworkListObjects!DllGetClassesAndCategories+000C38DC
mscorwks!CLRToCOMWorker+0000019A

Stack #2:

ntdll!RtlAllocateHeap+00000274
MSVCR90!malloc+00000079
MSVCR90!operator new+0000001F
DisplayCore!ClassFactory<CartographicLineSymbol>::CreateInstance+000000D8
AfCore!FactoryCache::CreateInstance+00000095
AfCore!FactoryCacheDispatcher::CreateInstance+00000050
AfCore!AfCreateInstance+00000040
SystemCore!ObjectStream::LoadObject+00000297
AfCore!LoadObject+0000005E
DisplayCore!MultiLayerLineSymbol::Load+000000F8
SystemCore!ObjectStream::LoadObject+000003E9
AfCore!LoadObject+0000005E
GraphicElements!LineElement::Load+00000076
SystemCore!ObjectStream::LoadObject+000003E9
AfCore!LoadObject+0000005E
GraphicElements!GroupElement::Load+00000082
SystemCore!ObjectStream::LoadObject+000003E9
AfCore!LoadObject+0000005E
AfCore!ExtraSafeReadObject+00000072
PageLayout!PageLayout::Load+00000174
SystemCore!ObjectStream::LoadObject+000003E9
MMStoredStuff!DllGetClassesForCategoryReg+00011B32
MMStoredStuff!???+00000000 : 2DB04186
MMStoredStuff!DllGetClassesForCategoryReg+0000D9E5
MMStoredStuff!DllGetClassesForCategoryReg+0000C297
mmPageTemplatesUI!DllCanUnloadNow+00020402
mmPageTemplatesUI!DllCanUnloadNow+00038EAF
mmFrameworkListObjects!DllGetClassesAndCategories+000C38DC
mscorwks!CLRToCOMWorker+0000019A

UPDATE #1:

Schneider Electric's support was able to reproduce the problem and has created a bug. I was told that if a fix was made at all, it would not be available until the 10.2 release.

I'll post more as I get more information.


texts: Calculating area difference between two polygon shapefiles using QGIS

texts:
I have two polygon shapefiles which contain about 30 features. These features have different boundaries in both shapefiles so the areas are different too.

I want to calculate that difference in the area for each feature using QGIS. Shall I use the intersection tool for that?


texts: Getting stripe across tiles without buffering image?

texts:
I am new to geotools.  I need to get a horizontal stripe of bytes from the image data without buffering the whole image.  I've looked at the various envelopes and pretty sure I could use just Rectangle2D since I don't need the coordinaterefsystem.  But for the life of me I can't figure out how to use the rectangle/envelope without doing:

GridCoverage2DReader gridReader = format.getReader(this.getInputFile(), hints);
GridCoverage2D gridCoverage = gridReader.read(null);

I don't see a replacement in JAI-EXT for ImageReadDescriptor.

For those wondering why, legacy code appears to have a memory leak. Below is the code where I see avail memory drop with each loop.  I used Thread.sleep(x) so I could watch it in top.  This is not a heap issue, the machine core dumps with an out of memory error and an exit code of 137.  Java mission control shows the heap is cleaned nicely.

RenderedImage img = ImageReadDescriptor.create(ImageInputStream iis, ...);
... //calculate stripeSize

byte[] writeBuff = new byte[stripeSize*numLines];
for(nCalculatedChunksOfLines){
  Rectangle tb = new Rectangle(0, n, sample, linesInChunk );
  Raster stripe = img.getData(tb);

  ...rest of processing bypassed for debugging....
  if(stripe != null){
    continue;
  }
}

I would rather replace the legacy code but if there are no simple options then a suggested fix to the legacy code would be helpful.  I'm going to test closing and opening the stream for each rectangle retrieve.


texts: Rotating labels in cartocss with controlling the rotation basepoint

texts:
I am trying to rotate multiple text labels in cartodb/cartocss. I am offseting each label so that they appear to be stacked using text-dx and text-dy. Using text-horizontal-alignment: left; and text-vertical-alignment: bottom; I can control the alignment for unrotated text, but when I use text-orientation: 90; (see documentation) it rotates around the middle and center of the string (i.e., text label).

So, when it is oriented at 90 degrees, it is no longer aligned left (or to the bottom of the page). Since my strings/labels are not all the same length, it causes unpredictable rotated label placement.

Any tips for how I can control the pivot/orientation point when the text is rotated/oriented? Or, is there a way that I can pivot around my marker?

The code I was playing with is here:

    #test_project_los::labels_sl
    [zoom >= 16]
    {
      text-name: [vol_sl];
      text-face-name: 'DejaVu Sans Book';
      text-size: 13;
      text-horizontal-alignment: right;
      text-vertical-alignment: middle;
      text-label-position-tolerance: 0;
      text-fill: #229A00;
      text-halo-fill: fadeout(white, 40%);
      text-halo-radius: 1;
      text-dy: -20;
      text-dx: -20;
      text-orientation: 90;
      text-allow-overlap: true;
      text-placement: point;
    }

    #test_project_los::labels_st
    [zoom >= 16]
    {
      text-name: 'Appears way right ' + [vol_st];
      text-face-name: 'DejaVu Sans Book';
      text-size: 13;
      text-horizontal-alignment: right;
      text-vertical-alignment: middle;
      text-label-position-tolerance: 0;
      text-fill: #3E7BB6;
      text-halo-fill: fadeout(white, 40%);
      text-halo-radius: 1;
      text-dy: -20;
      text-dx: -30;
      text-orientation: 90;
      text-allow-overlap: true;
      text-placement: point;
    }

    #test_project_los::labels_sr
    [zoom >= 16]
    {
      text-name: [vol_sr];
      text-face-name: 'DejaVu Sans Book';
      text-size: 13;
      text-horizontal-alignment: right;
      text-vertical-alignment: middle;
      text-label-position-tolerance: 0;
      text-fill: #FF2900;
      text-halo-fill: fadeout(white, 40%);
      text-halo-radius: 1;
      text-dy: -20;
      text-dx: -40;
      text-orientation: 90;
      text-allow-overlap: true;
      text-placement: point;
    }


texts: How can I remove easily vector leftovers within dissolved polygons?

texts:
I got an official government shapefile with zip code areas. After I dissolved the zip codes to the desired aggregated polygons, there are partial vector leftovers within the polygons. How can I clean them up systematically?

Below the before and after.

Before:

And here the dissolved map:

My version of QGIS is 2.18.1 with GRASS.


texts: Creating building footprints from LAS files

texts: I wish to create building footprints from the classification in our LAS files.  I have filtered out the building points, using the point cloud filter. But I now trying to join the outermost points of the buildings, which I can't figure out how to achieve this. I understand that this may be a little hit and miss with the LAS classifications but I am happy to accept this. I am a inexperienced FME user.


texts:
I have dGPS positions collected using a Leica 1200 system (base station and rover) recorded at the same site on multiple occasions over two years (being used to georeference terrestrial laser scanner data).

The positions were collected and post processed in WGS84 and converted to OSGB36 using OSTN15/OSGM15 as if they were ETRS89 using Grid InQuest.

As I understand it, OSTN15/OSGM15 should take ETRS89 lat/long as input and ETRS89 coordinates are slowly diverging from WGS84 coordinates.

Questions:
1) Will my positions be accurate relative to eachother, i.e. if I had surveyed the same mark each time, would the easting, northing and ODN elevation be the same or would it gradually shift due to tectonic plate movement? (I am interested in cm scale changes so the 2.5 cm annual change would affect me)

2) If not, can I correct them using a time dependent transformation?

3) Although I understand any transformation will introduce error, will any offset introduced by OSTN be consistent at a particular position?


texts:
There is a situation where I have annotation in dgn file as an attribute info, eg: it is circumference of trees trunk. I have a feature class in esri file geodatabase named ‘tree’ which has an attribute ‘tree circumference’. Is it possible to copy these annotation values  into the appropriate field of attribute table of feature class (point). Below is the picture of annotation layer and tree layer in dgn file. Here the ellipses and arcs are already converted into esri point feature class. Now the annotations associated with them are to be copied in attribute table. Is this possible with FME ? How ?

Result of match from NeighborFilter (max dist  = 5)

Possible solution is :

Use Sorter with 'min'. Next, to separate the first item(ie: minimum distance value) from the rest use DuplicateRemover. Since we dont need any values other than the nearest, use AttributeCreator to change them to blank.


texts:
I have a requirement where the user would like to see the output value in a raster as a decimal.  The user will input a series of ranges with a cost associated to them.  For instance

  0-5, 1.20

  5-10, 2.7

  10-15, 5.32

These ranges map to values on existing raster.  I want to create a new raster using the ReclassOp (or something else that I haven't found yet) that will map the ranges from the existing raster to the appropriate cost value.  For example, anything from the initial raster within the 5-10 range will be mapped to 2.7 in the new raster.  So my new raster ideally would have 3 values, 1.20, 2.7 and 5.32.  When the user clicks on an area in the raster in ArcMap, he should see one of these 3 values (1.20, 2.7 and 5.32) instead of the integer value (1, 3, 5).

This is my current code-

INumberRemap remap = new NumberRemapClass();
remap.MapRange(range.MinValue, range.MaxValue, (int)Math.Round(range.Value));
IReclassOp op = new RasterReclassOpClass();
output = op.ReclassByRemap(_source, remap, _mapNoData);

You can see that I have a value as a decimal (it is representing cost to the user), but INumberRemap.MapRange only accepts an int (currently I am rounding the double).  I searched the documentation and could not find a good alternative.  Is there a way to remap a range to a decimal value?

EDIT - Code with dividing the raster

My current problem is that this outputRaster (the return raster) has equivalent max and min values (equal to the min value).  When I look at the Statistics - Band_1 property in ArcMap on that raster I see "Statistics have not been calculated".

    //Create the remap with ranges
    INumberRemap remap = new NumberRemapClass();
    foreach (NumberRange range in Ranges)
    {
        int value = maintainDecimals ? (int)Math.Round(range.Value * 100) : (int)Math.Round(range.Value);

        if (range.IsNoData)
            remap.MapRangeToNoData(range.MinValue, range.MaxValue);
        else if (range.MinValue == range.MaxValue)
            remap.MapValue(range.MinValue, value);
        else
            remap.MapRange(range.MinValue, range.MaxValue, value);
    }
    IReclassOp reclassOp = new RasterReclassOpClass();

    // Reclassify the geodataset
    IGeoDataset reclassDataset= reclassOp.ReclassByRemap(_source, remap, false);

    //Create const geodataset
    RasterMakerOp rasterMaker = new RasterMakerOpClass();
    IGeoDataset constDataset = rasterMaker.MakeConstant(100.0, false);

    //Divide the reclassed raster to get decimals
    IMathOp mathOp = new RasterMathOpsClass();
    return mathOp.Divide(reclassDataset, constDataset);


texts: Correctly merging DEM rasters and eliminating edge effects

texts:
I'm using 10m DEM rasters from the NRCS Data Gateway for a large portion of California. The problem is some strange tiling that appears after merging, but isn't present in the original data (before merging). The grids create low points, so flow accumulation etc is drawn into them.

I've tried mosaic to new raster, create raster catalog and filter in ArcGIS 10.2, merge and build vrt in Qgis2.4 but the problem persists throughout.

The grids represent flat lows in the curvature raster on the merged edges of the DEM rasters. Note that those lines are the edges of BOTH of the mosaiced DEM, creating 2 lines.

The lines are flat lows, which creates a problem trying to do any sort of slope stability or hydrologic analysis...which is the whole point of this exercise.

Unmosaiced DEM Rasters:

Curvature raster generated from Mosaiced raster using bilinear interpolation:

EDIT: The data and mosaic method were the problems: The data itself has edge effects built in and I was using nearest neighbor instead of bilinear interpolation


texts:
I have a four band raster, and I would like to generate a new raster that has the argmax of the bands for each pixel. For example, if band 3 holds the maximum value at that pixel, the new raster should have a value 3 at that pixel.

I feel that this is a simple operation, but I am having a hard time with it. There is no argmax function in Postgres that I can find. And the callbackfunc syntax for ST_MapAlgebra() defeats me.

Is there a simple way to do this with an existing function?

Or are there examples of callback functions (along similar lines to this question, although that question wasn't really answered)?


texts: Clip layout view to output graphics extent using ArcPy?

texts:
I need to implement, within this very simple piece of code, the option that replicate the "Clip to output graphics extent" available from the ArcGIS 10.0 Export Map UI standard command. Is it possible? The purpose is to iterate a list of MXDs and export them their layout view to a graphic format.

mxd = arcpy.mapping.MapDocument(r"E:\XXX\\" + m)

arcpy.mapping.ExportToPNG(mxd, outDir + '\\' + str(m))

Please note I need to export the Layout View and NOT the Data View.

I have further investigated and I think the solution probably lays along the following lines:

extract from the DATAFRAME_ELEMENT its width and height
change the page printer size according to the values above using suggestions from geonet.esri.com/thread/59446
center the dataframe to the page margin
export the layout to a graphic file.

I will try the approach above and publish any results.


texts: How to exclude a certain attribute from being added to the Processing script output layer?

texts:
I have a working script to split a layer by field attributes and save the single features as files:

##[User scripts]=group
##input=vector
##class_field=field input
##output=output file

from qgis.core import *
from PyQt4.QtCore import *
from processing.core.VectorWriter import VectorWriter

layer = processing.getObject(input)
provider = layer.dataProvider()
fields = provider.fields()
writers = {}

class_field_index = layer.fieldNameIndex(class_field)

inFeat = QgsFeature()
outFeat = QgsFeature()
inGeom = QgsGeometry()
nElement = 0

feats = processing.features(layer)
nFeat = len(feats)
for inFeat in feats:
    progress.setPercentage(int(100 * nElement / nFeat))
    nElement += 1
    featAttr = inFeat.attributes()
    classField = featAttr[class_field_index]
    if classField not in writers:
        outputFile = output + '_' + classField + '.shp'
        writers[classField] = VectorWriter(outputFile, None, fields,
                                      provider.geometryType(), layer.crs())
    inGeom = inFeat.geometry()
    outFeat.setGeometry(inGeom)

    # here I'd like to exclude the 'classField' from the attributes in the output files
    outFeat.setAttributes(featAttr)

    writers[classField].addFeature(outFeat)

for writer in writers.values():
    del writer

However, I'd like to exclude the field that is used for the spliting of the source layer from the output files. The attributes are used for naming of the output files and thus aren't needed anymore..

Can anybody please provide some help for coding this?


texts: Preventing the "Coordinate Reference System Selector" query/window after processing.run()

texts:
I want to convert a DXF with python/ogr2ogr.
Unfortunately QGIS 2.99 always asks for the CRS.

QSettings().value('/Projections/defaultBehavior')
'useGlobal'
QSettings().value('/Projections/layerDefaultCrs')
'EPSG:25833'

processing.run('gdal:convertformat',{'INPUT':'d:/tar/mytest.dxf','OPTIONS':'','OUTPUT': 'c:/../x.shp'})

Without the processing-plugin the following code also produce a CRS Window

from processing.tools import dataobjects
parameters={'INPUT':'d:/tar/mytest.dxf','OPTIONS':'','OUTPUT': 'c:/../x.shp'}
context=dataobjects.createContext(None)
alg = QgsApplication.processingRegistry().createAlgorithmById('gdal:convertformat')
alg.checkParameterValues(parameters, context)

Is there a way to prevent the "Coordinate Reference System Selector" query/window?

I need this code for my plugin "Another DXF Importer / DXF2Shape Converter". This command is called 3-5 times per DXF. If there are 5 DXF-s, the dialog window appears 15-25 times. The query of the CRS is absolutely superfluous, since no transformation should take place and also no projection file must be written.

The code for QGIS 2.xx

processing.runalg('gdalogr:convertformat','d:/tar/mytest.dxf',0,None,'d:/tar/xx.shp')

works without Problems (without window)


texts:
As per the suggestions given by Michael to one of my earlier query (URL given below), I have tried to develop a sample code to perform rasterize using GDAL Java API methods.  I had written based on the Python sample from cookbook.

https://gis.stackexchange.com/questions/305164/using-gdal-proximity-function?noredirect=1#comment498883_305164

I am performing rasterize on a shapefile and the above code generates output but image in blank.

publicpublic void rasterizeShp() {

    //# Define pixel_size and NoData value of new raster
    int pixel_size = 25;
    int NoData_value = 255;

    gdal.AllRegister();
    org.gdal.ogr.ogr.RegisterAll();
    int bufferDistance = 20;

    String shpFilePath = "C:/img/GLXMABufferRadius.shp";

    //# Open the data source and read in the extent
    DataSource sourceDs = org.gdal.ogr.ogr.Open(shpFilePath,gdalconstConstants.GA_ReadOnly);
    Layer sourceLayer = sourceDs.GetLayer(0);

    SpatialReference sourceSrs = sourceLayer.GetSpatialRef();
    double[] extent = sourceLayer.GetExtent();

    //# Create the destination data source
    double x_res = ((extent[1] - extent[0]) / pixel_size);
    double y_res = ((extent[3] - extent[2]) / pixel_size);
    System.out.println("x_res -------" + x_res +  "y " + y_res );
    int xCor = (int)x_res;
    int yCor = (int)y_res;

    System.out.println("xCor -------" + xCor);
    System.out.println("yCor -------" + yCor);

    String output = "C:/img/myAttrOp2.tif";
    org.gdal.gdal.Dataset target_ds = gdal.GetDriverByName("GTiff").Create(output, xCor, yCor, 1, gdalconst.GDT_Byte);
    target_ds.SetGeoTransform(new double[]{extent[0], pixel_size, 0, extent[2], 0, -pixel_size});
    //target_ds.SetGeoTransform(sourceLayer.GetExtent());
    Band band = target_ds.GetRasterBand(1);
    band.SetNoDataValue(NoData_value);
    band.FlushCache();

    int[] intArr = {1};

    //# Rasterize
    gdal.RasterizeLayer(target_ds, intArr, sourceLayer, null);

}

Is this the right approach for rasterizing. Also I would like to know the correct way of setting geotransform.  I sometimes get "ERROR 1: Attempt to create 2x0 dataset is illegal, sizes must be larger than zero".

Please let me know if there are corrections to be made in the below code.


texts:
I work in an Oracle shop. My I.T. department uses ArcGIS Online to host data for use in data collection activities.

The data goes through a yearly cycle:

Export from our Oracle enterprise geodatabase
Import to ArcGIS Online
Edit on tablets or phones via the Collector app during the course of the year
Export to a file geodatabase
Process the data
Import back into the Oracle enterprise geodatabase

I'm involved in steps 1, 5, and 6. But, it has occurred to me that I have no idea what goes on in step 2. Out of curiosity, what is the underlying DBMS of ArcGIS Online?


texts: Using local mean sea level (LMSL) as altitude for converting latitude longitude altitude to ECEF

texts:
I am new in Geographic area and would like to convert latitude and longitude coordinates into ECEF x,y coordinates so that I could do some linear calculation of line or circle equation. I was able to get the elevation of the coordinate using Google Maps Elevation API and it says that the value is local mean sea level (LMSL). Can I use this value as altitude for ECEF? Or should I also use the geoid height value from EGM2008 or combine both LMSL+geoid height as altitude?

Case in point:

Latitude = 52.394033° N = 52° 23' 38.52" N
Longitude = 13.133462° E = 13° 8' 0.46" E
GPS ellipsoidal height (Google elevation) = 42.27206802368164 (meters)
Geoid height (EGM96) = 40.229 (meters)

Which one I should use as altitude for ECEF?


texts: Finding municipal-level administrative boundary data for India?

texts:
Does anyone know how to find a shapefile with boundaries of Indian municipalities?

I have tehsil boundaries as for year 2001, but I would like boundaries of municipalities and, ideally, also historic ones.

I believe Census of India has the information but I couldn't find any maps or shapefiles...


texts: There are some grid-based global addressing systems. NAC (Natural Area Code) is licensed (expensive), but I've seen one from Japan which used more interesting honeycomb-like divisions. I cannot find reference to the latter one, can anyone help with this?


texts: I am looking for a good resource (MOOC, book, etc..) to teach Python for remote sensing applications.  I am already familiar with http://www.rsgislib.org/ - are there any other resource out there for a beginning Python programmer?


texts: PyQGIS 3 - Add New Layer to Spatialite Database

texts:
For the most part, working with QGIS layers has been a pleasure, but I've run into some difficulties recently when trying to create a new vector layer in the SQLite database I've been using.

I've created a new database with the QgsVectorFileWriter like this:

writer = QgsVectorFileWriter(
    "/path/to/db/testdb.sqlite",
    "utf-8",
    myFields,
    QgsWkbTypes.PointZ,
    QgsCoordinateReferenceSystem(
        4326, QgsCoordinateReferenceSystem.EpsgCrsId),
    "SQLite",
    ["SPATIALITE=YES"])

And noticed that QGIS creates a default new table called "testdb" (named after the database) that works fine to load and save as to a layer.

I want to create a new layer with a new table now in the same database, so I went over to the docs at https://qgis.org/api/classQgsVectorFileWriter.html.

Using QgsVectorFileWriter()

I saw there was an ActionOnExistingFile parameter in the constructor to create a new layer & table, only to be disappointed when I noticed this parameter had been mysteriously dropped in the python implementation of the QgsVectorFileWriter, according to the PyQGIS docs.
Feeling hopeful, I threw in an action=QgsVectorFileWriter.CreateOrOverwriteLayer keyword, and even tried typing out all 11 arguments to avoid the keyword argument, but no luck.

# Doesn't work
writer = QgsVectorFileWriter(...
    action=QgsVectorFileWriter.CreateOrOverwriteLayer)

Using QgsVectorFileWriter.writeAsVectorFormat()

I did see that the python writeAsVectorFormat method had an option to specify the new table action, so I went ahead and wrote up something like this, using the SaveVectorOptions class:

options = QgsVectorFileWriter.SaveVectorOptions()
options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
options.layerName = "testlayer"
writer = QgsVectorFileWriter.writeAsVectorFormat(
    layer,
    "/path/to/db/testdb.sqlite",
    options)

But... The layer argument assumes I already have a valid layer from a valid data source. I need to create the source first before I can create the layer though, right?

Using a memory layer...?

If one can't create the data source first, I thought maybe loading a layer with a temporary "memory" provider and then writing it to the database would be an option. I'm not sure how I would specify the other parameters that I need for the vector layer (such as the fields and the CRS) using this:

layer = QgsVectorLayer("PointZ", "testlayer", "memory")

and then saving it to the database. I'm not even sure whether this method would work, but it seems a little hacky to me in any case. But I could be wrong...


texts: Is it possible to print one item legend item in italic in print composer?

texts:
Is it possible to have one item in italic font (species names) and the rest in "normal"?
I use QGIS 2.6.

The legend consists of two items which should be written normal (trap positions and vegetation) and one should be in italic, as it is convention that Latin species names are written in italic fonts.


texts: Adding CSV file with lat/long to map in different CRS in QGIS

texts:
I'm very new to QGIS, so I'm sure the solution to this is common sense.

I have a map in the WGS 84 zone 58S CRS and am trying to add lat/long coordinates from a text file to it. The issue: the WGS 84 58S CRS uses meters as its reference, but my coordinates are, like I said, lat/long. This is causing my points to show up in the wrong place. If I load the text file in a different CRS (WGS 84), they show up in the right place, but then the data is in a different CRS, which is no good.

How do I take lat/long data in a text file and import it into my maps CRS, and have the data points show up in the correct locations on the map?


texts: Adding label to QgsVectorLayer in PyQGIS?

texts:
My goal  is to define labels programmatically for the case similar for the figure below:

But I didn't find any guidance how to do so using Python.

I have a feeling it could be done using QgsLabel but didn't find any working exmaple how to do so.

I seek some example how to do so;

If I'm wrong with QgsLabel idea, can you explain briefly the right approach for this?

As usual some code example could be very useful.


texts: Fixing a corrupted geometry column with Python GeoPandas

texts:
I had a shapefile which I converted to a geodataframe, and created a geometry column. I did some manipulation with it in pandas. I am now trying to convert it back into a geodataframe and eventually a shapefile but cannot make it work. I have tried the usual methods:
import pandas as pd
import geopandas as gpd

df = pd.read_csv('DirectMail.csv')
df = gpd.GeoDataFrame(df, geometry='GEOMETRY__C')

and get the following error. This is the first row in the dataframe by the way:
TypeError: Input must be valid geometry objects: MULTIPOLYGON (((1241887.500623662 949342.9993114918, 1241827.375743911...
After looking at other Stack Overflow threads I tried this way using shapely wkt
import pandas as pd
import geopandas as gpd
from shapely import wkt

df = pd.read_csv('DirectMail.csv')
df['GEOMETRY__C'] = df['GEOMETRY__C'].apply(wkt.loads)

and I get this error:
WKTReadingError: Could not create geometry because of errors while reading input.
If I check to see if the geometries in the 'GEOMETRY__C' column are valid like this...
import geopandas as gpd
import shapely as sh

gdf = gpd.read_file("DirectMail.csv")

gdf['validation'] = gdf['GEOMETRY__C']

I have the following result:
0   False
1   False
2   False

and so on...
This makes me wonder, is the 'GEOMETRY__C' column in the wrong data type?  I checked the data type of this column and it is pandas dtype('O')
Does anyone have other suggestions what I can do?  Do I need to changes the datatype of the geometry column?


texts: Elevation obtained from USGS DEM is different from Lidar elevation

texts:
For a project I'm working on I need to calculate height of objects for which I know their spatial 2D coordinates (e.g. buildings). I'm using Lidar data to extract the highest Z value for an object, and then I subtract from it the ground elevation obtained from a DEM.

However, after getting some inconsistencies in heights calculated that way, I thought I should compare the ground level of the DEM I'm using with that of the Lidar data.

In all the tests I've done I found that the elevation I got from the DEM is larger than the ground level obtained from the Lidar data by about 2 meters.

I obtained both the DEM and the Lidar data from USGS. As the Lidar data points to use in the comparison I used points classified as ground (Classification = 2).

The area I used for the comparison is in Chicago, Illinois. The DEM I was using is for Cook county, it is in latitude/longitude coordinates. It is a raster tif file. The Lidar data is in EPSG:3435 CRS. Both are using, as far as I understand, GRS-80 ellipsoid.

I compared a sample of 1000 random points selected from the Lidar data, and I also compared points of a flat, open field that is about 200 x 200 square meters. I verified that the points I compared are of the same area. For the open field, the DEM data says its average height is 183.5 meters, the Lidar data says it is 181 meters. In both tests the height distribution was very narrow, with a standard deviation about 0.35.

The tests, and my project, are written in R.

I was expecting to get, more or less, the same value for ground elevation, and I don't understand why the values are so different. Am I missing something, or is something wrong with the method I'm using for the comparison?

Edit

I'm adding links to the specific Lidar and DEM data I was looking at:

Lidar
DEM


texts: Import OpenStreetMap tiles in a Geoserver Tiled Layer

texts:
I've got OpenStreetMap Tiles downloaded from OSM and stored in a repository, with the structure: Zoom/X/Y.png (own of OSM).
I need to import all of them in one of Geoserver directories in order to use them as a Geoserver tiled layer.

How can I create a layer like this?

Do I need some plugin?

Can I create a generic tiled layer and after replace all tiles with mine?


texts: Getting white lines in merged GeoTIFF files after ImageMosaicing on GeoServer?

texts:
I have merged 21 GeoTIFF images using GDAL utilities. Each image is the result of merging 500 smaller tiles. I am using ImageMosaic in Geoserver but I am getting a mosaic that contains white lines and distortions. The lines and distortions become smaller and finer when zooming in.

What could be the reason behind having those white lines?

The commands that I used are:

gdalbuildvrt -srcnodata 255 -vrtnodata 255 -a_srs EPSG:27700 -input_file_list tiff_list.txt mosaic.vrt

gdal_translate -of GTiff -co TILED=YES -co BIGTIFF=YES -co COMPRESS=JPEG -co JPEG_QUALITY=80 -co BLOCKXSIZE=512 -co BLOCKYSIZE=512 -co PHOTOMETRIC=YCBCR mosaic.vrt mosaic.tif

gdaladdo mosaic.tif -r average --config COMPRESS_OVERVIEW JPEG --config JPEG_QUALITY_OVERVIEW 60 --config INTERLEAVE_OVERVIEW PIXEL --config PHOTOMETRIC_OVERVIEW YCBCR 2 4 8 16 32 64 128 256 512


texts: Why does a TIF file lose projection information when a pixel value is changed?

texts:
I am trying to change the pixel value of a tif file (containing projection information) but after the pixel value is changed it is losing projection information.

I am using the Python imaging library for this and the code is like this

import Image
im = Image.open("some tif file")
im.putpixel((2,3),0)

So please help me in solving this ..


texts: CartoDB toggle buttons not working with my data?

texts:
We are taking CartoDB data for this Canadian election map -- http://cdb.io/1TabLyz -- to create a version that adds toggle buttons to filter by political party.

However, the buttons only work properly for data filtered out of the original map. Which would be fine, except the infowindows don't load in the new version.

Has anyone encountered anything similar? This is how it looks: http://www.cpac.ca/sites/leaders_map.html


texts: Aggregating detached polygons using open source software

texts:
Is there any open-source library to aggregate detached polygons?
I know that there is a tool in ArcInfo, but it is commercial.


texts:
I install Qgis web client on windows 7 and it's working fine!
Now I wont to setup search feature!

First I added line in Apache config file:
-LoadModule wsgi_module
modules/mod_wsgi.so
In GlobalOption.js I set
var useGeoNamesSearchBox = false;
var searchBoxQueryURL
="http://localhost:89/qgis-web-client-master/wsgi/search.wsgi?query=";
var searchBoxGetGeomURL
="http://localhost:89/qgis-web-client-master/wsgi/getSearchGeom.wsgi";

When i try to search, i do not get anything!
Is there any example for searching setup?


texts: I want to cluster houses into groups of 32 in a QGIS point file, with the clusters arranged to select houses which are in rows.  All my houses are in rows of 5-10 and I want to exclude houses which are not in rows from the clusters. Does anyone know of something I can use in QGIS?


texts: Allowing ArcGIS Server accessing our Geodatabase Enterprise over the web

texts:
In my organization there is an installation of ArcGIS Server that is linked with an Enterprise geodatabase (PostgreSQL) in a DMZ for editing only over the web.

However, we have another installation of ArcGIS Server that is linked with an Enterprise geodatabase (SQLServer) for internal use only. The main data of the organization are stored in the internal SQLServer, and its not accessible from the web.

We use some automatization tools to copy edited data from PostgreSQL to SQLServer and vice versa, but the desynchronization is becoming a problem.

The DMZ was designed to prevent putting our data at risk. However, we know there is many organizationd that don't really care about that, and are not afraid to expose their main database through the web.

We would like to study the possibility to let ArcGIS Server to have read and write access to our main database over the web and at the same time removing the DMZ.

What could be the issues of doing such change in our organization? Is there some way to attenuate the risks if such risks exist?


texts: Using deegree with Oracle DBMS without installing Oracle Spatial?

texts:
I'd like to use deegree for its WFS capacities connected to an Oracle 12 database where some geo information are available: stations location (latitude, longitude, altitude) + realtime information (parameters changing each 15 min).

In addition to licence costs for Oracle spatial, it seems to me useless to install it, since all necessary information can be extracted through SQL queries.

How can I use deegree with Oracle DB?


texts: How can I create 1 seamless tab file for master map?

texts:
This is my first post on here so please excuse if it's in the wrong place.

I've processed some mastermap files which display absolutely fine. However, during the processing it would not let me create just 1 seamless tab file as the index was over 2Gb. I therefore had to split it in to attributes, i.e. Topo_Point_Seamless, Carto_Area_Seamless etc.

Now although these work if you load them all separately, it makes it more awkward to work with. Ideally I'd like just 1 seamless tab file i can turn on and off in MapInfo and geoset's etc. However i've tried the seamless creator in MapInfo (11.5) but it says I cannot create a seamless file from other Seamless files.

Does anyone know how i could go about this?

Many thanks,

Nick


texts:
I am working with Arcmap 10.1 and census tract level data for 1990 and 2011. I am trying to compare different options for tracking changes between the tracts.

The tracts change between the census polls, and there are some tracts from 1990 that merged to form a single tract in 2011, and those that split up.

I would like to know what are some options I have to keep track of the changes in the tracts and how I can compare demographic changes from 1990 to 2011?


texts: how to show wfs and wms sources in boundlessgeo gxp viewer

texts: I'm using boundlessgeo gxp viewer. I have a lot of different wms and wfs link and want to show those remote sources in my map viewer application. I have a geoserver and basic layers in it. I can add my local sources normally. What is the way of adding external wms and wfs sources?


texts: Error While Converting UTM data to Lat-Long Coordinate

texts:
I have a a dataset of Vancouver,BC where Coordinate values are projected in UTM Zone 10 . Now  I can convert it using spTransform by the following piece of code and looks like the it correctly transforms the coordinates ,as I can see if I plot them in a map .

cord_transform <- function(df){
  coordinates(df) <- ~X+Y
  proj4string(df) <- "+proj=utm +zone=10 +ellps=GRS80 +units=m +no_defs"
  df.lat.long <- spTransform(df,CRSobj=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
  df.lat.long

}

However  I realized I dont need a SPdataframe ,rather my purpose serves if I can transform the UTM coordinates in lat -long coordinates in the original data.table (since the data is huge I want to use data.table ) but when I follow  this answer,and write a similar code for the data.table

func <- function(x,y){
  proj4string <- "+proj=utm +zone=10 +ellps=GRS80 +units= m + no_defs"
  inp <- list("x"=x,y="y")
  print(inp)
  pj <- project(inp,proj4string,inverse = TRUE)
  print(pj)
  latlon <- list(lat = pj$y,lng = pj$x)
  latlon

}
crime.data.2016.jan[,c("lat","lng"):=func(X,Y)]

I get this error :

Error in project(inp, proj4string, inverse = TRUE) :
  unknown unit conversion id

a subset of the data looks like this :

         NEIGHBOURHOOD        X       Y
 1:       Renfrew-Collingwood 497445.9 5454283
 2: Central Business District 491580.1 5458745
 3:                           492757.5 5458792
 4:              South Cambie 491371.1 5455550
 5:                  West End 490417.0 5459169
 6:               Shaughnessy 489926.6 5454845
 7:                   Marpole 490409.0 5450869
 8:            Mount Pleasant 492545.4 5456743
 9:                  Fairview 490767.8 5456820
10:                 Kitsilano 488624.5 5457182
11:                 Killarney 497326.3 5452640
12:          Hastings-Sunrise 496763.8 5458220
13:  Kensington-Cedar Cottage 494917.4 5455246
14:                Kerrisdale 488711.0 5452722
15:       Victoria-Fraserview 495207.1 5452375
16:                Strathcona 493345.6 5458515
17:         Dunbar-Southlands 486589.5 5454941
18:                  Oakridge 490997.4 5453021
19:                    Sunset 492937.0 5451823
20:             Arbutus Ridge 488145.0 5455054
21:        Grandview-Woodland 494940.8 5458213
22:                Riley Park 492607.3 5455102
23:           West Point Grey 485649.1 5456957
24:              Stanley Park 490531.0 5460906
25:                  Musqueam 485508.0 5453015
                NEIGHBOURHOOD        X       Y


texts: How do I merge layers using an external Python script to call the processing module?

texts:
This is what I'm trying to run:

import sys, glob

# Prepare the environment
from qgis.core import * # qgis.core must be imported before PyQt4.QtGui!!!
from PyQt4.QtGui import *
app = QApplication([])
QgsApplication.setPrefixPath("/usr/bin/qgis", True) # The True value is important
QgsApplication.initQgis()

# Prepare processing framework
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
driver2 = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open("file1.shp", 1) #1 is read/write
dataSource2 = driver.Open("file2.shp", 1) #1 is read/write
layer1 = dataSource.GetLayer()
layer2 = dataSource.GetLayer()

Processing.runAlgorithm("qgis:mergevectorlayers", layer1, layer2, '/home/name/Desktop/output.shp')

But I get this error:

[Finished in 0.2s with exit code 1]
[shell_cmd: python -u "/home/none/Desktop/testattribute/merged.py"]
[dir: /home/name/Desktop/testattribute]
[path: /home/name/torch/install/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games]

I'm using Processing 2.12.2.

Now I'm running this code from @gene with the addition that I'm trying to use it in a loop that iterates through a list of shapefiles:

meta = fiona.open(filelist[0]).meta
    for shapefile in filelist:
        with fiona.open(exportname, 'w', **meta) as output:
            for features in fiona.open(shapefile):
                output.write(features)

But I do not end up with a valid shapefile.  What am I doing wrong?

I ended up sorting it out into a somewhat hacked looking script.

def mergeshp(firstshp, secondshp, index):
    exportname = 'mergeshp_' + str(index) + '.shp'
    meta = fiona.open(firstshp).meta
    with fiona.open(exportname , 'w', **meta) as output:
       for features in fiona.open(firstshp):
           output.write(features)
       for features in fiona.open(secondshp):
           output.write(features)
    return exportname

And I called it with this code:

numfiles = len(shapelist)

i=0
mergedshp = mergeshp(shapelist[i],shapelist[i+1], i)
i += 2
while i<numfiles:
    mergedshp = mergeshp(mergedshp, shapelist[i], i)
    i += 1


texts: How to display URL / email in popup window

texts:
I would like to display data from Geoserver in URL format using popup window. I created a popup window that displays the data as string. What would be the best way to display email as a link to outlook and website as url?

    var url = standorte
                .getSource()
                .getGetFeatureInfoUrl(
                    evt.coordinate,
                    map.getView().getResolution(),
                    map.getView().getProjection(),
                    {
                        'INFO_FORMAT': 'application/json',
                        'propertyName': 'firma,adresse,telefon,email,web'
                    }
                );

    reqwest({
        url: url,
        type: 'json',
    }).then(function (data) {
        var feature = data.features[0];
        var props = feature.properties;
        var info = "
" + props.firma + "

" + props.adresse + "
" + props.telefon + "" + props.email + "<c>" + props.web + "</c>";
        popup.show(evt.coordinate, info);
    });

}

});


texts: Distance matrix returns points in wrong order QGIS 3.4

texts:
I'm using QGIS 3.4 on Windows 10 and scripting in pyqgis.  I have a layer representing building centroids that is in a projected CRS. I want to find the azimuth between each building and its nearest neighbour in the same layer, so the input and target points being in the correct order is important.  However, when I run the following script the closest pairs are output each time with the points in the same order rather than reversed.  This renders the output useless.  I have tried changing the parameters for MATRIX_TYPE but have kept the other parameters as shown.  The data is a vector point layer in projected format.

import os
import ogr

from qgis.core import (
    QgsVectorLayer,
    QgsProject,
    QgsProperty
)
import processing

QgsProject.instance().removeAllMapLayers() # remove all of the layers from the active map

Building_centroids = QgsVectorLayer('C:/Users/ADMIN/Desktop/Mod04AssessedEx/data/bldg_centroids.shp', "Building_centroids", "ogr") # create vector layer, ogr (or something else) needed

if not Building_centroids.isValid():
    print("Layer failed to load!")

QgsProject.instance().addMapLayer(Building_centroids) # add layer to map

distance_matrix = processing.run("qgis:distancematrix",

{
'INPUT': Building_centroids,
'INPUT_FIELD': 'OBJECTID',
'TARGET': Building_centroids,
'TARGET_FIELD': 'OBJECTID',
'MATRIX_TYPE': 0,
'NEAREST_POINTS': 1,
'OUTPUT': 'memory:'
}
)['OUTPUT']

QgsProject.instance().addMapLayer(distance_matrix) # adds the results as a temporary layer to the active map

# this calculates the azimuth between each structure and its nearest neighbour

for structure in distance_matrix.getFeatures():

    structure_geometry = structure.geometry()
    print (structure_geometry)
    structure_geometry = structure.geometry().asMultiPoint()# this line unpacks the multipoint geometry into an iterable or indexable list

    print(structure_geometry[0])
    print(structure_geometry[1])

    azimuth_2_neighbour = structure_geometry[0].azimuth(structure_geometry[1]) # is this working?? no buildings are due west of their nearest neighbour, or is that coincidence

    print(azimuth_2_neighbour)


texts:
Is it possible to edit POSTGIS tables directly in QGIS 2.6.1.   I've been using The DB MANAGER but can only really do DDL type edits....Or create/import a brand new layer.   I did see a SQL Window but that would be too complex for my customers to navigate.  I believe in the past I've seen a Plugin or something to do this but from what I've read the DB Manager is the plugin of the future for Postgis.

An alternative is to pull into QGIS, Save as a shapefile and then re-publish to POSTGIS.  This approach doesn't work well in our environment.

I'd like to be able to edit PostGIS directly via some type of Table Edit (similar to what I can do with shapefiles in ArcGIS).