texts: Overlaying Google Map in standalone PyQgis application?

texts:
I have created a standalone PyQGIS app (within Python itself, not a plugin) with PyQt and PyQGIS API. It generates my map really well. I would like to overlay a Google Map layer as base map just like what can be done within QGIS Desktop itself. Is this possible? If so, how should I go about doing it? Should I use QWebView or something?

Following is a snippet of my code:

#  Importing necessary classes and libraries like PyQt4 (for making GUI), qgis for Python-Qgis connection
#  Some classes are useful for overlaying OpenLayers.
#  Refer https://github.com/joystor/BaseMap2Image/blob/master/baseMap2Image.py
import sys
import signal
from PyQt4 import QtGui
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.gui import *
from qgis.core import *
from qgis.utils import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
# QtWebKit is used for application that interacts with the web
from PyQt4.QtWebKit import *
from functools import partial

# Specify App
app = QtGui.QApplication([])
# Specify which folder is QGIS installed
QgsApplication.setPrefixPath("/usr", True)
# Opening QGIS
QgsApplication.initQgis()
# Show setting of parameters to check if it matches that on QGIS (log messages at the bottom right -> general tab)
#print('QGIS App Settings')
#print QgsApplication.showSettings()
# Call main window of Qt GUI
main_win = QtGui.QMainWindow()
# Set up Qt GUI frame
frame = QtGui.QFrame(main_win)
main_win.setCentralWidget(frame)
# Specify Qt GUI layout (grid)
grid_layout = QtGui.QGridLayout(frame)
# Call up a shapefile to a vector layer
# Load vector layer
data_source = "/home/victorzhiyulee/Desktop/PJU Only Yeah/PJU Only Yeah Split 2_2.shp"
layer_name = "dun"
provider_name = "ogr"
layer = QgsVectorLayer(data_source, layer_name, provider_name)
# Quick check if layer is valid
if layer.isValid():
    print('Layer is valid')
else:
    print('Layer is invalid')
# Display window
main_win.show()
#lst = layer.id()
# 0.5 - Defining QGis Map Canvas
map_canvas = QgsMapCanvas()
grid_layout.addWidget(map_canvas)
map_canvas.setCanvasColor(QtGui.QColor(255, 255, 255))
#list_of_layers.append(lst)
# Add the big layer into registry
QgsMapLayerRegistry.instance().addMapLayer(layer)
print('Count of Layers: ', QgsMapLayerRegistry.instance().count())
print('Layers included are')
for lapis in QgsMapLayerRegistry.instance().mapLayers():
    print(lapis)
Senarai_layer = QgsMapLayerRegistry.instance().mapLayers().values()
print 'Or in QGIS language: ', Senarai_layer
for lapisan in Senarai_layer:
    print(lapisan.id())

# Setting up map canvas
canvas_layer = QgsMapCanvasLayer(layer)
map_canvas.setLayerSet([canvas_layer])
map_canvas.zoomToFullExtent()

# Printing
img = QImage(QSize(800, 600), QImage.Format_ARGB32_Premultiplied)
color = QColor(255, 255, 255)
img.fill(color.rgb())

p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)
render = QgsMapRenderer()
render.setLayerSet([layer.id()])

rect = QgsRectangle(render.fullExtent())
rect.scale(1.1)
render.setExtent(rect)

render.setOutputSize(img.size(), img.logicalDpiX())

render.render(p)
p.end()

# Show the window
main_win.show()
# Running the script
app.exec_()


texts: Knowing if map has specific custom control in OpenLayers 3?

texts:
I created my own custom control for openlayers 3 as follows:

window.app = window.app || {};
var app = window.app;

/**
 * Control Creado para moverse por el mapa (quitar controles de dibujar...)
 */
app.Move = function(opt_options) {

  var options = opt_options || {};

  var button = document.createElement('button');
  button.innerHTML = '<i class="fa fa-hand-o-up"></i>';

  var this_ = this;

  function move_ (){
      console.log();
      map.getControls().forEach(function (control) {
          console.log(control, 'control', window.app.GetFeatureInfo(), 'gfi');
          if(control instanceof (window.app.GetFeatureInfo())) {
              control.activar(false);
          }
          else if(control instanceof (new app.Draw())){
              control.removeDraw();
              control.removeModify();
              control.activar(false);
          }
      });
  }

  button.addEventListener('click', move_, false);

  var element = document.createElement('div');
  element.setAttribute('data-toggle', 'left');
  element.setAttribute('title', 'Mover');
  element.setAttribute('data-content', 'Desplazarse por el mapa');
  element.className = 'move ol-unselectable ol-control';
  element.appendChild(button);

  ol.control.Control.call(this, {
    element: element,
    target: options.target
  });

};
ol.inherits(app.Move, ol.control.Control);
map.addControl(new app.Move());

In this control I want to 'turn off' other custom constrols.
What I have done is to get all controls in the map and look up for my specific control using instanceof method.

  map.getControls().forEach(function (control) {
      //console.log(control, 'control', window.app.GetFeatureInfo(), 'gfi');
      if(control instanceof (window.app.GetFeatureInfo())) {
          control.activar(false);
      }
      else if(control instanceof (new app.Draw())){
          control.removeDraw();
          control.removeModify();
          control.activar(false);
      }
  });

But i'm getting this error:

TypeError: invalid 'instanceof' operand window.app.GetFeatureInfo(...)

How can I get a valid instance of custom control to compare it with the same custom control setted in the map?


texts: Adding shapes to polygon feature in the same row

texts:
I've got somehow simple issue but cannot process it.
I am working with shapefile in QGIS (3.10 on MacOS), where having several polygon features stored as rows. I want to somehow (manually) edit some of those features (pink one on the print screen) to add additional one (pink one the print screen) or two polygons nearby. To have them together in the same row, so just reshape polygon somehow as I need to process with area calculation and other analysis and need them to be in the same row.
Tried to google but maybe I am missing something obvious


texts: Using Gazetteers to mine locations from text

texts:
I am starting to work with Gazetteers and trying to find something that will actually "mine" geospatial locations from sets of text. Academically, this is a difficult task as texts are based on context, and as such, things like THE BRIDGE, would be very difficult to identify. The Bridge is the Golden Gate, FYI(See how my second sentence lets you know where I am).

The end result that I am looking for is the ability to feed text into a tool and output locations with XY's. (Potentially with some level of "confidence", which would be asking for a lot).

I am trying to build a tool to feed on twitter and produce geospatial outputs(To visualize political opinions for academic purposes).

I have some experience with Python.
I know this is a very complex question looking for a very simple answer.But any pushes in the right direction would be very helpful.  And although this seems like a text-mining/data mining question, I believe it is also a GIS question, because it is surrounded by the concept of context, the idea of assigning some absolute space to a relational thing, as well as there is a matter of spatial accuracy, and precision related to the idea of assigning absolute coordinates to things.


texts: How to download Planet 3-month mosaics using the API

texts:
I would like to download Planet 3-month mosiacs using the RestFUL API, but I'm not sure what to specify for the item_type. The documentation here https://developers.planet.com/docs/api/items-assets/#assets lists a bunch of item types but they seem to be individual tiles.

Is it possible to use Planet's API to download 3 month mosaics, and if so, how?


texts: How to implement search tool get values in dropdown ArcGIS JavaScript 3.x

texts:
How to get point?

When search button clicked  i am trying search button click to go  to point in ArcGIS JavaScript my code is this.

How to push the values into search button?


texts: Adding result of geoprocessing service -- consume in ArcMap desktop

texts:
I am new to ArcGIS Server.  I have a Python script tool that I've published as a geoprocessing service.

The script takes a point feature class, writes a subset of the input FC as a new FC in the scratchworkspace, creates some new fields and calculates values for these fields based on user input parameters.

The service and all of the calculations the script is performing, works fine - I can go into the scratch workspace on the server, find the unique job ID folders, and see my results.

Eventually, I want to consume it via the Flex API, but for now, I am trying to get it working by consuming it in ArcMap Desktop.

My problem is getting these results loaded into the map.  I've followed the advice given in response to this question.  I've added an output parameter,  tried to load the results as a FeatureSet and set the output parameter to be the FeatureSet.

When I do this using the following code, I get an error that I "cannot load table into a FeatureSet".  This doesn't add up because my output is a point feature class (I have confirmed this by looking at it in the scratchworkspace on the server).

result = "%scratchworkspace%\Scratch.gdb\output"
res = "in_memory/result"
r = arcpy.CreateObject("FeatureSet")
r.load(res)
arcpy.SetParameter(9, r)

Any advice or direction would be appreciated.  Thanks,
-Erik


texts:
I was trying to add a label and a primitive over the same area like so:

var viewer = new Cesium.Viewer('cesiumContainer');
var scene = viewer.scene;

var instance = new Cesium.GeometryInstance({
  geometry : new Cesium.RectangleGeometry({
    rectangle : Cesium.Rectangle.fromDegrees(-130.0, 40.0, -115.0, 60.0),
    vertexFormat : Cesium.PerInstanceColorAppearance.VERTEX_FORMAT
  }),
  attributes : {
    color : new Cesium.ColorGeometryInstanceAttribute(0.0, 0.0, 1.0, 0.8)
  }
});

viewer.entities.add({
        position : Cesium.Cartesian3.fromDegrees( -122.19, 46.1914),
        label : {
            text : 'Some Text',
            font : '24px Helvetica',
            fillColor : Cesium.Color.WHITE,
            outlineColor : Cesium.Color.BLACK,
            outlineWidth : 2,
            style : Cesium.LabelStyle.FILL_AND_OUTLINE
        }
    });

scene.primitives.add(new Cesium.Primitive({
  geometryInstances : [instance],
  appearance : new Cesium.PerInstanceColorAppearance()
}));

The end result - the text appears under the image like so:

How can I add Label and primitive so that they don't hide one another?
Am I missing something here?


texts: How to write a custom canvas item for SVG marker by subclassing QgsMapCanvasItem in QGIS?

texts:
I want to subclass QgsMapCanvasItem for SVG marker in QGIS, I looked up a similar topic in stackexchange and QGIS mailing list, I found that this is similar to mine, I did exactly what @Nathan W did. But when I instantiated custom canvas item, the app broke down to a break point. It says:

  "This may be due to a corruption of the heap..", see more in error.jpg.

.

I don't know why, can anyone help me out of it?

Here is my code:

//main.cpp
int main(int argc, char ** argv)
{
QgsApplication app(argc, argv, true);

QgsMapCanvas * mypMapCanvas = new QgsMapCanvas();
mypMapCanvas->enableAntiAliasing(true);
mypMapCanvas->setCanvasColor(Qt::white);
mypMapCanvas->freeze(false);
mypMapCanvas->setVisible(true);
mypMapCanvas->refresh();

QString fileName = "D:\\qgis\\svg\\symbol\\poi_tower_communications.svg";
TFMapCanvasSvgItem mTFSvgItem(mypMapCanvas, fileName);
mTFSvgItem.setPos(0, 0);

return app.exec();
}

//TFMapCanvasSvgItem.h --custom canvas item named TFMapCanvasSvgItem
class TFMapCanvasSvgItem: public QgsMapCanvasItem
{
public:
    TFMapCanvasSvgItem(QgsMapCanvas* canvas, QString fileName = "");
    ~TFMapCanvasSvgItem();

    void paint(QPainter* painter);
    QRectF boundingRect() const;
    void setFilePath(const QString& file);
    QString filePath() const {return mFilePath;}

private:
    QSvgRenderer mSvgRenderer;
    QString mFilePath;
};

//TFMapCanvasSvgItem.cpp
TFMapCanvasSvgItem::TFMapCanvasSvgItem(QgsMapCanvas* canvas, QString fileName): QgsMapCanvasItem(canvas), mSvgRenderer(fileName)
{
}

TFMapCanvasSvgItem::~TFMapCanvasSvgItem()
{
}

void TFMapCanvasSvgItem::paint(QPainter* painter)
{
if (!painter)
{
    return;
}

QRectF viewBox = mSvgRenderer.viewBoxF();
if (viewBox.isValid())
{
    setPos(this->toCanvasCoordinates(QgsPoint(0.0, 0.0)));
    mSvgRenderer.render(painter, viewBox);
}
}

QRectF TFMapCanvasSvgItem::boundingRect() const
{
QRectF viewBox = mSvgRenderer.viewBoxF();
qreal halfWeight = viewBox.width() / 2;
qreal halfHeight = viewBox.height() / 2;
return QRectF(-halfWeight, -halfHeight, viewBox.width(), viewBox.height());
}

void TFMapCanvasSvgItem::setFilePath(const QString& file)
{
mFilePath = file;
mSvgRenderer.load(mFilePath);
}


texts: How to view more than one layers balloon description box in Google Earth

texts:
This may be a lot simpler than I am making it, but I have racked my brain over this for some time:

I have two shapefiles in ArcMap -
The first is a county wide layer made up of hundreds of postcode sector polygons, each showing an average house price. This has been categorized and colour coded into seven categories.

The second is a layer with a select few regions within the county (as polygons) which are thought to be good locations for a new development which I want with no fill but a black outline.

Now I have exported these using the layer to kml tool with no problems and they keep the symbology and formatting I want, however when looking at the description bubble when I click on an area it only shows me the information for the average house price data.

Is there a way that I can view both that and the name of the location thought to be a good for a new development, whilst keeping the formatting of both?


texts:
I have a large amount of data in Excel spreadsheets.

I have been converting that data into Attribute tables to use in ArcMap.

The data includes:

names of locations
geo-location
descriptions of each location

I need to find an easy way that users, with only basic knowledge of ArcMap, can query the information in the tables and plot only the locations that meet the query parameters on the map, excluding the remaining data.

What would be a good way of going about this, without having to write a script?


texts: Polygons not displayed properly. Grey colour above the data

texts:
I am trying to import a shape file cotaining Wind data over Sweden. I have to zoom in very much to be able to view the data under the grey colour that appears. I figure it has to do something with the border of the polygons but after having had them disabled the problem still remains.

If I zoom out with the globe button I first get the desired format (See picture: https://i.imgur.com/qbqkcOO.jpg) but as the whole loading is completed the gray layer appears above it (See picture: https://i.imgur.com/316guZP.jpg) and it is making it impossible for me to begin analysing the data.
GIF of how it appears when I zoom out: https://i.imgur.com/2Xj5iRB.gifv
Why does it appear grey and how can I view it as portrayed in the first photo?

I am rather new to GIS so this could be a very "easy" fix.


texts: Order matched geometries on line by distance in PostGIS

texts:
I want to be able to find Points near (with a given threshold) a given Line-Geometry and Polygons that it crosses.
I already have that part of the problem solved, but now I'd like to order the results. Fist match first, last match last.
Here's a sample:

My result should be something like:
--+-------------------------
1 | point A
2 | polygon A
3 | point B
4 | point A

I've tried ST_Line_Locate_Point which seems to do basically that, but this will only give me the closest match, right?
So my problem would be the following. My input line crosses a certain point twice. One time a bit closer on the way back e.g. a bit further away.
Is there anything build into PostGIS to help me with that? Or will I have to do the analysis of something like that on my own?
How would you approach this problem?


texts: Correcting class codes of LiDAR files in ArcGIS

texts:
Is there any procedure to automatically correct lidar classification? My LAS Dataset has water listed as ground and buildings listed as vegetation.

Just wondering if there is any option to fix this.


texts: How to publish Qgis2threejs output online?

texts: I want to have something like this: http://maps.navigasi.net/examples/stack_layer/soil_layer.html , my problem is because when i put the Html file outside of the original directory it does not run.


texts: Finding density centroids within polygons in PostGIS

texts:
I want to find the highest density of OpenStreetMap nodes within a dataset of census geometries. I want to tie the census data to the location of the node which is most central in terms of population. This is required for a drive time analysis. I am aware that street nodes are not completely representative of population but it is the best we have and I need to run drive times.

I have a dataset of census geometries and assigned the id of the polygon geometries to the nodes that are covered by each polygon. I use st_coveredby to assign the blkgrp (US Census block groups) id to each OSM node.

I would like to run a query which takes all nodes within a 1000m buffer into account and calculates the density for that node based on the relative distance to all the other nodes in the 1000m buffer. The query will only look at nodes that are within the same block group.

So far I get the count using st_dwithin as well as an aggregate for the distances (either sum or avg)...

Now I just need to work out how to aggregate the distance to get the best density value for the points. I can then find the point with highest density and move the census area centroid to that location.

UPDATE osm_nodes_oh
  SET c_grp = Q.c
  FROM (
    SELECT
      t1.id i,
      count(t2.id) c,
      sum(t1.geom <#> t2.geom) d_sum,
      avg(t1.geom <#> t2.geom) d_avg
    FROM osm_nodes_oh t1, osm_nodes_oh t2
    WHERE st_dwithin(t1.geom, t2.geom, 1000)
      AND t1.blkgrp = t2.blkgrp
    GROUP BY t1.id) Q
  WHERE osm_nodes_oh.id = Q.i;


texts:
I am having trouble getting .shp layers to "line up" on QGIS.

Now, I opened them both in ArcGIS and noted there that they had the same projection. I opened them up in a text editor, and just to be sure, copied the data from one .prj file into the other. Reopened them in Arc and had no problems whatsoever.

However in QGIS I can't get them to come into the same map. I've tried every combination of CRS / Layer menu options that I can think of, nothing seems to work.

They seem to have the same Projected Coordinate System, but different Geographic Coordinate Systems: this I assume from going to the full extent, and note that both layers are on the map but hundreds of KM apart.

I tried setting them to the same GCS in QGIS properties, but then they both vanished. ("Zoom to Layer didn't work, they were gone.)

Is there are process I should be following here?

If the data were corrupted, wouldn't duplicating the .prj data be enough to fix the problem? If it's a GCS issue, how do I change that?

Mostly I'm looking for a process I can apply every time I import data to QGIS so I'm guaranteed to work around any CRS issues.


texts: Obtaining the Feature History from OpenStreetMap

texts:
I've seen research papers (Haklay 2010) using the data regarding the number of edits (or the contributors working on) on the features of OSM. I am using QGIS and download the area I want to study. However, data obtained this way are limited to shapefiles. How could I systematically obtain the number of edits made on each feature of OSM? This information have been used as a proxy for the data quality of OSM, which is something I want to evaluate.

References:

Haklay, 2010, How Many Volunteers Does it Take to Map an Area Well? The Validity of Linus’ Law to Volunteered Geographic Information


texts:
Forester, not a geospatial or IT professional, currently using QGIS and SAGA for generally county level geoprocessing of multipspectral imagery and DEM's (cluster analysis, zonal statistics, stream channel modeling).

Currently exceeding practical processing limits for laptop/windows configuration using HP EliteBook w/4GB Ram, Intel Core i5, 32 bit Win7 OS. Need advice on logical next step in hardware/OS/software.

Future ambitions include deployment of PostGIS and derivation of canopy height models from stereo matched imagery. Thinking low end workstation and Ubuntu OS.

I'm a one man shop trying to operate a bit beyond my capabilities. I do not speak fluent IT or GIS, need a balance of horsepower, idiot proof, moderate learning curve and intermediate support. Cost is important but functionality, flexibility and capability are primary.


texts: Use Python to find points based on 'month' and 'year' fields

texts:
I have a shapefile based on 27 locations, with historical values at each location for each month of the year going back to 1974. The file has 11,000 values. I have been trying to work on a script to select the 27 locations for each month of each year to eventually create shapefile for each monthly set of the 27 locations.

import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C://project"
LocationFile = "C://project//27location.shp"

FieldList = arcpy.ListFields("27location.shp")

for field in FieldList:
    print field.name
    # check for fields

cursor = arcpy.da.SearchCursor(Locationfile, ["NAME", "VALUE", "MONTH", "YEAR"])
for row in cursor:
    # select based on monthly values at the 27 locations into a seperate shapefile.

I hope this makes sense I have been struggling to interpret code.


texts: Map frame in the map composer doesn't show the project layers when the size of the map is too large. Small size maps work fine

texts:
I drew a map in QGIS print composer, but it does not render my layer from the QGIS project. When I reduce the size of my map it works. Here are two pictures to illustrate my issue.

Why is that and how can I change it?


texts:
I'm using IDLE. GDAL command works in QGIS' Python console but not Python 2.7 (ArcMap) and 3.8 (separate install) IDLEs.

So far, I have tried

Downloading this (GDAL/OGR IN PYTHON) and following the instructions in the README text file. Did not work.
Looked at the GDAL releases here (Stable Releases) but there is no GDAL release for Python 3.8.1 MSCv 1916. I followed this article (Installing GDAL for Windows) for this step.
Installed GDAL through the OSGEO4W installer but when I ran a code with this line, it did not work for Python 2.7 and 3.8 (but worked in QGIS python console)

from osgeo import ogr


texts: Drawing "proportional" buffer using QGIS or GRASS?

texts:
I'm working on aquatic vegetation and have a shapefile representing floating vegetation in a lake. By some studies, we know that submerged vegetation extends for a certain extension based on floating vegetation (e.g. submerged vegetation extends 3 times floating vegetation).
I'd like to draw potential submerged vegetation extension, so I'd need something like a "proportional buffer" drawing a polygon slight where the original polygon is slighter and larger where the original polygon is larger, and the resulting buffer distance has to be the result of the width of the original polygon multiplied by a given value (e.g. 3).
I know this is difficult, but does exist a way of doing it?
I attach an image where the brown polygons represent the floating vegetation and the red polygon, drawn by eye, should represent the result I'm looking for (more or less)


texts: Recommended way to update standalone QGIS on Windows?

texts: QGIS Standalone installer suggests a default installation location like c:\Program Files\QGIS Wien, including version name in the path. Does this mean the recommended way to update to a new version would be to uninstall the old version, and then install a later version? Or can we just as well install a newer version on top of the previous version? (in which case it makes more sense to change the location to e.g. c:\Program Files\QGIS).


texts: ITRF2008 to WGS84 Datum Transformation in ArcMap?

texts:
If I have data in IRTF2008 as shapefile is it possible to change this to WGS84 and how?

The need for this comes when trying to display the data in a third party program that doesn't recognize the .prj for the IRTF 2008 data.

I'm Using ArcGIS 10.2.2 and can't seem to find a transformation for this application.


texts: Role of boundaries when using shapely within function

texts:
According to the shapely documentation:

  object.within(other)

  Returns True if the object’s boundary and interior intersect only with the
  interior of the other (not its boundary or exterior)

from shapely.geometry import Polygon

p1 = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
p2 = Polygon([(0, 0), (1, 0), (1, 1)])
print(p2.within(p1))
# True

But p2 boundary intersect with p1 boundary ((0,0), (1,0)).


texts: .NET console application running super slow all of the sudden

texts:
I have a standalone .NET console application that (through ArcObjects) processes DWG files, extract coordinates, and pushes them to SQL Server. It's smoking fast and can process a DWG in a matter of seconds. It's been running fine in production for quite a while now (over a year). All of the sudden last Wednesday it started running dead slow (40x longer). My server had a major power failure on the Monday night before this all started going wrong. When the server came back up, iSCSI drives were missing, all kinds of stuff was jacked up and it look 4-5 guys in infrastructure hours to get it all back up. I just found out about my process running slow today.

My question is this: Does anyone have any idea what I could look for as the culprit in something like this. The process runs, just really, really sloooooow. Like its CPU starved. If it was permissions, it would just fail. I don't expect to get an answer from my infrastructure team. Any insight is greatly appreciated.

EDIT: We just realized that ArcGIS Server hasn't written to any log files since this hard crash. The plot thickens...

Another edit: The drive this executable resides on is connected through iSCSI. I don't know a whole lot about iSCSI, but we are wondering if this is a communications issue between the server blade and iSCSI (which I think communicates via a fiber card). I have moved the .NET application over to one of the physical drives on the blade and we'll see what happens...


texts: PopupTemplate Content strips HTML Data-* Attributes in ArcGIS API for JavaScript 4.16

texts:
Working with the ArcGIS JavaScript API version 4.13 and was able to build popup templates that included data attributes to work with the Semantic UI library for their Tab module (https://semantic-ui.com/modules/tab.html). I recently upgraded to version 4.16 and now the template is broken. Upon examining the html when rendered it strips key data-attributes for tabs (data-tab) prevent it from working. Is there a way to force explicit content and prevent the template from being modified in version 4.16 or if this expected going forward?
Example (Click on the map to launch the popup)

Version 4.13 Working: https://codepen.io/booshwa/pen/ExPJgRz
Version 4.16 Broken: https://codepen.io/booshwa/pen/qBbwaJB


texts: How can I generate irregular grid containing minimum n points?

texts:
Given a large (~ 1 million) sample of unevenly distributed points - is it possible to generate irregular grid (in size, but could also be irregular in shape if that is possible?) that will contain specified minimum amount of n points?

It is of less importance for me if genereted 'cells' of such grid contain exactly n number of points or at least n points.

I'm aware of solutions like genvecgrid in ArcGIS or Create Grid Layer in QGIS/mmgis however they will all create regular grids which will result in an output with empty cells (smaller problem - I could simply discard them) or cells with count of points less than n (bigger problem since I'd need a solution to aggregate those cells, probably using some tools from here?).

I've been googling around to no avail and am open to both commercial (ArcGIS & extensions) or free (Python, PostGIS, R) solutions.


texts: How to ensure the values set in the CSV are captured within the GeoJSON feature values when merging CSV attributes into a GeoJSON

texts:
I create a GeoJSON but the merged values are null. The length of the datasets do not match -- is this a potential issue? If so, how do I merge for all covid matches. (FYI this isnt Covid19 related COV is short for covenant- housing).
I have tried the method suggested here. https://gis.stackexchange.com/a/359055/177082
import pandas
import geopandas
import fiona

gdf = geopandas.read_file(r'.\\UL1.geojson') # geojson file
df = pandas.read_csv(r'D:\ATMS\data\Blank.csv') #CSVfile

hvista=gdf.merge(df, on='Covid', how='left')

hvista.to_file("Hvista.geojson", driver="GeoJSON")

print(hvista)

and assuming it is a data type issue i have tried the following:
    import pandas
import geopandas
import shapely
from geopandas import GeoDataFrame

gdf = geopandas.read_file(r'.\\UL1.geojson') # geojson file
df = pandas.read_csv(r'.\data\Blank.csv') #CSVfile

gdfa=gdf.merge(df.astype(str), on='Covid', how='left')
gdfa=GeoDataFrame(gdfa, crs="urn:ogc:def:crs:EPSG::3857")
gdfa.to_file("Hvista.geojson", driver="GeoJSON")

print(gdfa)


texts:
I'm trying to calculate if a canon will be able to shoot over a height in the terrain.

lets saytThe canon have to be 15'000-18'000m from the target, and the impact degrees have to be over 30 degrees.

So the position of the target will be known, and the result i'm looking for are the areas where one can make the shot from. It would be fairly easy if one knew both the position of the target and the canon, but my trouble is that the canon is mobile and i do't know any efficient methods to calculate for all the possible positions...

I thought of making a spheroid or something and see how/where it would intersect with the elevation model, but i don't have the competence for it.

Tried to visualize me problem and desired output.

1 image: dilemma if the canon is to close to reach over the mountain

2 image: idea of result, red areas are areas where the shot can not be made from, while the green areas are ok

I'm using ArcGis for Desktop 10.3, also have 10.2 available, all extensions.


texts: count band numbers from multiband raster image

texts: My raster image has 10 bands (10 bands are 10 year), cell values of each one are 1, 2, 3. One is negative anomaly, 2 is no change area, 3 is positive anomaly. Now I want to find how many years out of 10 have the negative anomaly,how many years have positive anomaly etc. For example, the cell(10;10) has a value of 3 for all 10 bands, then the value of output raster in cell(10;10) should be 10. if cell(11;11) has a value 3 for 7 bands, then output value would be 7 etc.


texts:
I want to use the package movecost on R (first-time user). I have created a DTM and imported a point as a starting point, and a layer of points representing destinations. They all have the same CRS.
When I run it, I get this message:
Error in .local(.Object, ...) :

I ran the traceback() tool, and get the message down below, which I can't interpret (as I said, first-time user).
13: .local(.Object, ...)
12: initialize(value, ...)
11: initialize(value, ...)
10: new("GDALReadOnlyDataset", filename, silent = silent, allowedDrivers = allowedDrivers,
        options = options)
9: rgdal::GDAL.open(object@file@name, silent = TRUE)
8: .readRasterLayerValues(x, 1, x@nrows)
7: .local(x, ...)
6: getValues(mask)
5: getValues(mask)
4: .local(x, mask, ...)
3: raster::mask(accum_final, dtm)
2: raster::mask(accum_final, dtm)
1: movecost(dtm = areaDTM, origin = Start, destin = Destinations, move = 8,
       breaks = 0.5)

Here is my code:
library(lattice)
library(movecost)
library(raster)
library(rasterVis)
library(rgdal)
library(ggplot2)

setwd("/Users/Desktop/Data File/ANALYSIS")

# Load destination points from shapefile

Destinations <- shapefile("/Users/Desktop/Data File/ANALYSIS/Sites/Destinations.shp")

# Load starting point from shapefile

Start <- shapefile("/Users/Desktop/Data File/ANALYSIS/Starting_points/Start1.shp")

# Load raster DTM

areaDTM <- raster("/Users/Desktop/Data File/ANALYSIS/DEM/DTM.tif")

# plot on map

plot(areaDTM)
  plot(Start, add=TRUE)
  plot(Destinations, add=TRUE)

# use movecost package

result <- movecost(dtm = areaDTM, origin = Start,
                   destin = Destinations, move = 8,
                   breaks = 0.5)


texts: QGIS: JSON (GeoJSON) from Public API is empty if done by protocol

texts:
This Question is actually related to this one:
QGIS: JSON Load error: Invallid Layer cannot be loaded in: Origin is GeoJSON
Because the context is the same but I am going on this by a different angle.
Nevertheless let me repeat myself here:

To give some context I am testing a system, where we can get a JSON
more specifically a GeoJSON from a API. For now we are testing it with
the public API. Link to the Public API documentation is here:
http://api.gipod.vlaanderen.be/#!index.md

Now I am trying to retrieve this GeoJSON file by Portocol and the selecting GEOJSON as type of protocol in QGIS 3.4.5.
The problem is that the GeoJSON I get clearly is a table without geometry and its also empty.
This is the resource URL:
http://api.gipod.vlaanderen.be/ws/v1/workassignment/
Clearly there are records there, so why is this one empty?


texts: Shapefile of Baden-Württemberg and Rhineland-Palatinate with district boundaries

texts: Can anyone help me find shapefiles of German states Baden-Württemberg and Rhineland-Palatinate (including their districts)?


texts: Is it possible to define the field displayed on the Identify window in ArcMap?

texts:
I have a series of layers that are linked through relationship classes.  When the user clicks on a watershed (HUC = Hyrdologic Unit Code) with the identify tool a list of species and impaired rivers/streams is brought up.

The problem is the field that ArcMap decided to show for the impaired river/stream is not helpful at all in this context (some permanent ID).  I'd like it to show the actual impairment as defined by the Clean Water Act 303(d) dataset, which is in this same table.
Is there a way to define the field that ArcMap displays in the identify box??


texts: QGIS: Create line vector between line without drawing

texts:
I have transected the Transect line(white line) from a densified pink line
Now I want to draw the green line between my transected line.
Instead of doing it manually, is there another way to do it by using some processes of copy the line and offset?
TQ


texts: The interpolation of the Z value when snapping in the creation of a new feature is done well in QGIS 3.14, but when editing elements this does not happen. Is there any way to do it?


texts: How to do a Loop with a condition in Model Builder?

texts:
Can I loop a set of actions in model builder while a specific features rows are 1 or more.

Some background:

I have a point layer where I've created a workflow that:

first selects the row with the minimum value in a field containing distances to a line.

by using: [DISTANCE] in (SELECT min( [DISTANCE] ) FROM PointLayer) as my selection query.

I then create a buffer around this one selected point.

I plan to copy and append this point into some sort of result layer.

I then select by location and delete all points from the original PointLayer that are within the buffer with 'delete rows' or 'delete features' tool.

Now, since the original PointLayer is modified and some points have been deleted I want to restart the process and loop it until no points/rows are left in the PointLayer.

A new row will be the minimum since the old one was deleted in the previous run.

So the core question again, how can I make such a while loop in model builder and how do I apply the condition while rows > 0, I found there is a get count tool that counts rows.


texts: Can QGIS manage non-spatial PostgreSQL database?

texts:
I have a database with a table with plots. This table has a spatial field for the area. Also, I have another table with the owners. An owner can own some plots, and there is an "id_owner" field in the plot table.

My question is if it is possible to use the table "owner" in QGIS? If I can show the data, or make queries or update fields in "owner" table.

I have searched before but I have not seen exactly this question.


texts: Python Toolbox Value Table Parameter - getting separate outputs for different column inputs

texts:
I'm using a Value Table as a parameter for a Python toolbox so the user can input two related attributes. I'm struggling with how to use the values the user inputs as separate values/parameters. The code for my value table parameter looks like this:

    target_elements = arcpy.Parameter(
        displayName = "Target Element(s) Name",
        name = "element_name",
        datatype = "GPValueTable",
        parameterType = "Required",
        multiValue = True,
        direction = "Input")
    target_elements.columns = [['GPString','Target Element'],['GPString','EO ID (if known)']]
    target_elements.filters[0].type = 'ValueList'
    target_elements.filters[0].list = sorted(list(et_all.values()))

This code is giving me the desired toolbox operation. However, if the user input is "Abagrotis anchocelioides" for the first column and "53910" for the second column, the resulting output is a string that puts them together like this:

'Abagrotis anchocelioides' 53910

I'm not sure how to separate these values so that I can use them separately. I can't split on spaces because the first value often contains two words.

Is there another way to do this, or should I be using a different parameter type?


texts: QPyNullVariant not recognized outside of QGIS

texts:
I'm getting an error about an operation involving a QPyNullVariant:

TypeError: unsupported operand type(s) for *: 'QPyNullVariant' and 'float'

I wanted to put a check into my code to skip these values. However, my IDE (Eclipse) [and more over the terminal Python interface] doesn't recognize this as part of the PyQt4.QtCore module whereas it seems fine within QGIS. I have updated all my packages and restarted my computer.

import PyQt4.QtCore as pyqt
import sys
print(sys.version)    #1
print(pyqt.__file__)  #2
pyqt.QPyNullVariant   #3

in all cases this produces the following out put for 1 and 2:

2.7.6 (default, Jun 22 2015, 17:58:13)
[GCC 4.8.2]

/usr/lib/python2.7/dist-packages/PyQt4/QtCore.so

However in Eclipse (with PyDev) and in the terminal interpreter the call to PyQt module for a QPyNullVariant attribute fails:

AttributeError: 'module' object has no attribute 'QPyNullVariant'

whereas in the QGIS python interpreter you get:

<class 'PyQt4.QtCore.QPyNullVariant'>

Thoughts?


texts: Seeking free tile server like mapbox studio?

texts: Is there any free tile service/server like mapbox studio so that I can upload my .shp files and it will work like it does in mapbox studio tiles and get geojson tiles on Leaflet.js map?


texts: Is it possible to cache the same WMS service based on a vendor parameter?

texts:
I'm exploring optimization techniques for caching data coming back from a Digital Globe service. We have two unique identifiers we use for their WMS GetMap service that will give us low or hi-res raster images. I'd like to cache both, but the only different between the two is simply a vendor parameter called 'connectid'.

In the end, we'll cache at least one of them, but we'll need to make a choice of which one if it's not possible to do both.


texts: Lambert conformal conic isn't actually conformal?

texts:
Mathematically a conformal map has the property that the Jacobian of the transformation is a scaled version of a rotation matrix.  The transformation involved in a spherical Lambert conformal conic map projection doesn't have that property.  In what sense is the LCC conformal?

Details

Map projection in pseudocode (x,y are the map coordinates); I'm using the notation as described on pg. 104 of Map Projections - a Working Manual

x=rho*sin(theta)
y=rho0-rho*cos(theta)

Jacobian (phi,lam refer to latitude, longitude respectively):

J =   (drho/dphi)*sin(theta)     rho*cos(theta)*(dtheta/dlam)
     -(drho/dhpi)*cos(theta)     rho*sin(theta)*(dtheta/dlam)

This would be a scaled rotation matrix if (drho/dphi)==rho*(dtheta/dlam) however, since

 rho=R*F*cot( phi/2+pi/4)**n
 theta=n*(lam-lam0)

(R,F,n, lam0 are constants derived from the parameters defining the projection.) we have

drho/dphi=-n*R*F*cot( phi/2+pi/4)**(n-1)*csc( phi/2+pi/4)*0.5
         =-n*rho*csc(phi/2+pi/4)**2*tan(phi/2+pi/4)*0.5

rho*(dtheta/dlam)=n*rho

Thus, the Jacobian is not a scaled version of a rotation matrix, therefore the LCC is not a conformal map.


texts:
I want to union a set of Public Land Survey System (PLSS) polygons in a PostgreSQL/PostGIS layer, based on two fields, ‘sec_qtr’ and ‘qq’. A sample subset of polygons looks like this:

I want to select the polygons that have an identical ‘sec_qtr’ value and ‘qq’ values of ‘NE’ or ‘SE’ (i.e., the ‘East Half of the Northeast Quarter of Section 27, T1S, R1W’ in PLSS parlance). The resultant unioned geometry that I want is this:

However, the PostgreSQL script that I’m using:

SELECT st_union(geom) FROM work.tmp_qq1 where sec_qtr IN
(SELECT sec_qtr
  FROM work.tmp_qq1
  GROUP BY sec_qtr
  HAVING COUNT(DISTINCT(sec_qtr)) > 1)
AND qq = 'NE' OR qq = 'SE';

gives me this result:

What do I need to change in my SQL script to get both ‘01S01W27NE’ polygons unioned together?


texts: Can't get correct number of populated places for a place

texts:
I am trying to get a number of populated places for example for Andorra since it is a small number.  Here is a geonames site with 10 places:

http://www.geonames.org/AD/largest-cities-in-andorra.html

But if I do this query:

http://api.geonames.org/search?featureCode=PPL&username=genadinik&country=AD&style=full&maxRows=100

Only 4 out of the results have population.  Am I doing something wrong?  Is there a more accurate way of getting populated cities there?

Thanks!


texts:
I have GPS locations with date. These are points.

I want to convert the points to lines, sorted by date.

Using PostgreSQL

Tried:

SELECT observations.id,
       ST_MakeLine(observations.geom ORDER BY observations.date) AS newgeom
FROM   observations
GROUP BY
       observations.id
;

The problem is that the single line values gives the same coordinates. So, there will be no line. 0 and 1 are the same.

How is this possible?


texts: Non-default icon of OpenLayers Vector Feature

texts:
thats propably some minor problem, but i have deadline above my head and i can't get any clue whats wrong with my code, so I decided to come here for help. Im changing from markers to vector layer and i can't make my portal use any sort of non-default icon, whatever i put i externalGraphic style attribute doesnt have effect on map. I just see annoying orange circles.

Thats my code for layer and style:

var StylIkony = new OpenLayers.Style({
externalGraphic : '${symbol}',
graphicWidth    : 15,
graphicHeight   : 15,
cursor        : 'pointer'
});

var StylWarstwyIkon = new OpenLayers.StyleMap ({
default: StylIkony
});

warstwa_ikon = new OpenLayers.Layer.Vector("Ikony Lokali");
warstwa_ikon.StyleMap = StylWarstwyIkon;
map.addLayer(warstwa_ikon);

And this is already executed code from creating point features:

WspolrzedneIkony = new OpenLayers.Geometry.Point(2279231, 7127620);
Ikona = new OpenLayers.Feature.Vector(WspolrzedneIkony,
    { "symbol": "grafiki/IkonyLokali/10.png", "idLokalu": 1 });
warstwa_ikon.addFeatures([Ikona]);

WspolrzedneIkony = new OpenLayers.Geometry.Point(2279245, 7127630);
Ikona = new OpenLayers.Feature.Vector(WspolrzedneIkony,
    { "symbol": "grafiki/IkonyLokali/6.png", "idLokalu": 2 });
warstwa_ikon.addFeatures([Ikona]);


texts: Intersection of a MultipolygonZ and a floating point raster

texts:
I am trying to isolate the pixels in a floating point raster (loaded in PostGIS) that intersect my MultipolygonZ table (also loaded in PostGIS) and whose values are less than or equal to the z values of the MultipolygonZ table - essentially a 3d intersection.

I am running PostGIS 2.1.5. SELECT postgis_full_version() returns:

"POSTGIS="2.1.5 r13152" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER"

Here's what I have tried so far:

1) ST_Intersects(rast, geom) - this simply calculates a 2d intersection between the polygons and the raster - the z values of the MultipolygonZ are ignored.

SELECT (ST_DumpAsPolygons(rast)).*
FROM pgz_tab, dem
WHERE ST_Intersects(rast, geom)

2) ST_3dIntersects(rast, geom) - invalid arguments - apparently this only works with points and linestrings - http://postgis.net/docs/manual-2.1/ST_3DIntersects.html

SELECT (ST_DumpAsPolygons(rast)).*
FROM pgz_tab, dem
WHERE ST_3dIntersects(rast, geom)

3) Converting the MultipolygonZ to a floating point raster and then subtracting the two rasters. I can't seem to get the arguments right with the ST_AsRaster function - http://postgis.net/docs/manual-2.1/RT_ST_AsRaster.html even after reading A guide to the rasterization of vector coverages in PostGIS.


texts: How do I convert a .bil (and associated .clr, .blw, .hdr, and .stx) files into georeferenced rasters in R? OR to geotiffs using gdal_translate?

texts:
The title says it all.  I've got a bunch of .bil files, which I can open in R using the raster command pointed to the path to the .bil file.  But I don't know how to efficiently incorporate the contents of the other files, which seem to be metadata, into the raster layer in order to georeference it.

Alternatively, could I convert these sets of files to geotiffs using a system() call to gdal_translate?  What would be the syntax for that?


texts:
Does anybody know where I can find a database containing all (or a large percentage of) building addresses in the U.S.? I'm particularly interested in commercial properties.

Perhaps there's a google api I could query for this information?

I was able to find this question but it's about finding a database full of building footprints whereas here all I'm concerned about is finding a database full of addresses. It is also five years old.


texts: GeoServer 2.9.2 database connection pooling/thread locking issue

texts:
I am mocking GetMap requests. Scenario is that 50 users use only "database layers". Each user "sends" 8 GetMap requests every ~22s for ~300s which sums up to 5200 GetMap requests for vector data. Problem is that with these parameters I see that thread count starts to grow almost immediately and usually reaches ~2k threads. When inspecting with VisualVM I saw that 1.5k threads where BLOCKED: ~50% by org.geoserver.catalog.ResourcePool.getCacheableFeatureType(ResourcePool.java:943) and other half by org.geotools.jdbc.JDBCDataStore.getPrimaryKey(JDBCDataStore.java:1062). Thread count tends to grow when Java heap reaches max point but not always - it also grows rapidly even when heap is at ~50%.

When inspecting Tomcat logs I see:

org.geoserver.platform.ServiceException: Rendering process failed
...
Caused by: java.lang.RuntimeException: Unable to obtain connection: Cannot get a connection, pool error Timeout waiting for idle object
...

Increasing max connections from 20 to 100 in GeoServer makes threads accumulate even faster.

There are 6 CPU's and 5GB RAM 2 of which are given to GeoServer's Tomcat.

Any ideas on what is the root of such problem?

SOLUTION: DB server had not enough CPU power so sessions were waiting to be processed. Increasing CPU count from 2 to 8 seem to solve extremely high tread count issue.


texts: Arcgis 10.0 desktop: how to clear All transformations from map

texts:
I have a custom extension which adds some custom transformations to the map. When turned off, the extension is supposed to remove all the added transformations.
However, when calling IGeoTransformationOperationSet.RemoveAll(), nothing happens. The transformations are removed from the set but nothing changes in the map - the transformed features stay in the same place and do not move back to their original position.

I'm using Arcgis desktop 10.0 and my code is written in c#.

How can I turn off/on the custom transformations? What am I missing?
Thanks.


texts:
I have been doing this project for the last few weeks on QGIS3.4 with no problems and now every time I try and go on it QGIS does not respond and it does not load. I clicked on another project which loaded fine and this error appeared:

  An error has occurred while executing Python code:  RuntimeError:
  wrapped C/C++ object of type ContourGeneratorAlgorithm has been
  deleted

Does this have something to do with my project is now not loading?

Not sure if it is to do with a plugin bug.


texts: import SIP module. QGis 2.12 win 10 64bit. By book: Building mapping applications with QGIS

texts:
i have bought the book "building mapping applications with qgis" from eric westra.

Currently i'm following the steps described in the book. But at the verry first trial of an external app i get an error.

Traceback (most recent call last):
  File "C:\OSGeo4w64\apps\qgis\python\qgis\__init__.py", line 26, in <module>
    import sip
ImportError: No module named sip

I found several sip folders in my "osgeo4w64" folder.

further more i used the advanced osgeo installer?

as variables i setup:

osgeo4w_root = C:\OSGeo4w64
path = C:\Python27\;C:\Python27\Scripts;C:\ProgramData\Oracle\Java\javapath;C:\Program Files (x86)\Intel\iCLS Client\;C:\Program Files\Intel\iCLS Client\;%SystemRoot%\system32;%SystemRoot%;%SystemRoot%\System32\Wbem;%SYSTEMROOT%\System32\WindowsPowerShell\v1.0\;C:\Program Files (x86)\Windows Live\Shared;C:\Program Files\Intel\Intel(R) Management Engine Components\DAL;C:\Program Files\Intel\Intel(R) Management Engine Components\IPT;C:\Program Files (x86)\Intel\Intel(R) Management Engine Components\DAL;C:\Program Files (x86)\Intel\Intel(R) Management Engine Components\IPT;C:\Program Files (x86)\Intel\OpenCL SDK\2.0\bin\x86;C:\Program Files (x86)\Intel\OpenCL SDK\2.0\bin\x64;C:\Program Files\Condusiv Technologies\ExpressCache\;C:\Program Files (x86)\Microsoft SQL Server\100\Tools\Binn\;C:\Program Files (x86)\Microsoft SQL Server\100\DTS\Binn\;C:\Program Files (x86)\PFPS\system;C:\Program Files\SIMDIS\bin\amd64-nt;C:\Program Files (x86)\Boundless\OpenGeo\jetty\jre\bin;C:\Program Files (x86)\Boundless\OpenGeo\bin;C:\Program Files (x86)\Boundless\OpenGeo\sdk\bin;C:\Program Files (x86)\Boundless\OpenGeo\pgsql\bin;C:\Program Files (x86)\Skype\Phone\;C:\Program Files\FME\;c:\OSGeo4w64\apps\qgis\bin
pythonpath = C:\OSGeo4W64\apps\qgis\python
qgis_prefix = c:\OSGeo4w64\apps\qgis

as you can see the whole list at path has already multiple path's.


texts:
I am using leaflet.
I have a number of circles, and I want to imply a global event listener link for those circles which will add some other circles in one click to any circle, and remove those added circles by the second click.
How is that possible?

//the circles which will have an event listener
var circle1 = L.circle([44.6663888888889, 6.775],{radius: 150, color:'white',weight:.5, opacity:1,fillColor: '#ff9900',fillOpacity:1});
var circle2 = L.circle([44.6502777777778, 6.79333333333333], {radius: 150, color:'white',weight:.5, opacity:1,fillColor: '#008000',fillOpacity:1});
...

//create an onClick event targeting circle 1 and 2 to add/remove circle3 and circle4
map.on('click', function(e){...?...});

var circle3 = L.circle([44.6191666666667, 6.81222222222222], {radius: 150, color:'white',weight:.5, opacity:1,fillColor: '#008000',fillOpacity:1});
var circle4 = L.circle([44.6761111111111, 6.795], {radius: 150, color:'white',weight:.5, opacity:1,fillColor: '#008000',fillOpacity:1});
...


texts: Issue when transforming an EE.Image into numpy array

texts:
I use Google Earth Engine with Python 3 with Colab's notebooks.

I would like to transform the Google Images into Numpy arrays to be used for further processing.

I am using the toList Reducer to have a list of values corresponding to B8 band of a Sentinel2 tiles. This is a one dimension list so I need to reshape it properly.
However, I have a problem when doing a reshape on a large tile. Indeed, the output image contains a visual shift when the scale is too low (100m). When I do my tests with a much larger scale (1000m), this offset no longer exists.
Obviously I would like to have this result with much smaller scales.

There is the ouput for the 100m processing (KO) :

And there is for the 1000m (OK) :

You can find an example of the code here to test it for yourself (with scale of 1000 and scale of 100) :  https://colab.research.google.com/drive/1YFWxy3JOWiJ9UAaWsCVRvOwxx_0mobsd

import ee
import numpy as np

ee.Authenticate()
ee.Initialize()
 ############### Function ######################
 #---- Compute the transformation Earth Engine to numpy
def prepareData(collection, testingScale):
    CRS=collection.first().projection().getInfo()['crs'] # Get CRS
    img=collection.median().float()  # Compute an image and cast it into float
    countRes=img.reduceRegion(reducer=ee.Reducer.count(),\ # Get the size of the image
                 crs=CRS,\
                 geometry=collection.geometry().bounds(),\
                 maxPixels=1e13,\
                 scale=testingScale)
    img = img.reduceRegion(reducer=ee.Reducer.toList(),\ # Transform earth engine Image into a one dimension array
                 crs=CRS,\
                 geometry=collection.geometry().bounds(),\
                 maxPixels=1e13,\
                 scale=testingScale)
    dataB8 = np.array((ee.Array(img.get("B8")).getInfo()))
    dim=np.sqrt(countRes.get('B8').getInfo()).astype(np.int16) # Get the dimension for the reshape into two dimension array
    nrows, ncols = dim, dim
    image = dataB8.reshape((ncols,nrows))
    return image

############### Main body ######################
#---- Get an example tile of Sentinel2
point=ee.Geometry.Point([-70.066310,-21.341900])
collection = ee.ImageCollection("COPERNICUS/S2").filterBounds(point)\
                                  .filterDate("2018-01-01","2019-01-10")\
                                      .filterMetadata("CLOUDY_PIXEL_PERCENTAGE","less_than",10)\
                                  .select(['B8'])
#---- Prepare a numpy image
image = prepareData(collection, 1000) # Reprojection with scale 1000 OK !     Scale 100 non ok !
display(image)

#----- Plot the numpy image
import matplotlib.pyplot as plt
plt.imshow(image)
plt.show()

Has anyone ever had this kind of problem or had any ideas for explanation or improvement?


texts: How to determine ring orientation from lat/long coordinate points?

texts: We have some set of point coordinates in lat/long which contribute to a polygon. From this set how do we determine the polygon is clockwise or anticlockwise?


texts: Google Earth Engine: Supervised Classification fusion table error

texts:
I am currently doing a supervised classification (just two classes for detection of large-scale solar plants) in GEE. When I am using training polygons drew and create as polygons everything works just fine but when I use the same training polygons as a fusion table (made from KML) I keep getting this message. Here is the sample of the code.

var trainingData = solarDenmark;

// Train Sample Data
// Select Bands
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B11', 'B12', 'nd', 'NDWI'];
var input = sentinelComposite.select(bands);
var classifierTraining = input.select(bands)
      .sampleRegions({
        collection: trainingData,
        properties: ['solar'],
        scale: 30
      });

// Instantiate a CART and train it:
var classifier = ee.Classifier.randomForest(10).train({
  features: classifierTraining,
  classProperty: 'solar',
  inputProperties: bands
});

//Classify the image
var classified = input.select(bands).classify(classifier).clip(region2);


texts: When carrying out spatial point pattern analysis on my point files I get an error for standard distance

texts:
How to fix an error with standard distance in spatial point pattern analysis. When performing spatial point pattern analysis QGIS 3.4 is only carrying out the mean point calculation correctly and there is an error for the standard distance calculation. I am using EPSG 29902 (Irish National Grid) How do I get this to work.

See below for error.


texts: ESRI JS API multiple groupByFieldsForStatistics columns and having clause

texts:
I'm trying to query some data using ESRI JS API version 4.0.
The data table that looks like a longer version of this:

+-------------------------------------+
| Object_ID | Item_ty_co | Serv_grade |
+-----------+------------+------------+
|     1     |    "MH"    |     5      |
|     2     |    "MH"    |     3      |
|     3     |    "IN"    |     3      |
+-------------------------------------+

I am trying to group this data based on the item_ty_co and Serv_grade fields, and count the number of records in each group.

In regular SQL, my query looks like:

SELECT Item_ty_co, Serv_grade ,count(Object_ID) FROM public."pointAssets"
GROUP BY Item_ty_co ,Serv_grade
HAVING ITEM_TY_CO not in ('PN','GN','RN','CN')

which gives me exactly the result I want. However, I want to recreate this in ESRI JS API by querying the layerview like this:

const query = layerView.layer.createQuery();
query.outStatistics = [
   {
        onStatisticField: "Object_ID",
        outStatisticFieldName: "numberOfAssets",
        statisticType: "count"
   },
];
query.groupByFieldsForStatistics = ["Item_ty_co ","Serv_grade" ];
query.having = "ITEM_TY_CO NOT IN ('PN','GN','RN','CN')"

However, I'm having two problems here.

It seems like the API just ignores any columns in the groupByFieldsForStatistics array after the first element, and doesn't use them for the groupby. When I switch the items around, it just uses the first (0th) element.
query.having does not filter out the groups that I specify.

So the output for the above is an array of objects like this, with one record for the unique ITEM_TY_CO's:

0: Object { numberOfAssets: 4176, ITEM_TY_CO: "PN" }
1: Object { numberOfAssets: 4662, ITEM_TY_CO: "MH" }
2: Object { numberOfAssets: 5472, ITEM_TY_CO: "CN" }

I don't know where my mistakes are. I can write workarounds, but I'd really like a clean solution.

Querying the layer in the ArcGIS REST Services Directory interface with exactly the same parameters works perfectly. So I'm not sure why it doesn't work in the JS API. Also, the layer "supportsHavingClause" is set to true.


texts: Maximum Likelihood Classification in ENVI

texts:
I am a new ENVI user and I am trying to perform a supervised classification of a Landsat 8 image with the maximum likelihood algorithm.

In order to run the classification, I have to insert a Scale Factor Unit that, for uncalibrated integer data, is equal to the maximum resolution of the sensor (that, for a 16-bit instrument, is 65535).

In this case it works fine, but if I try to perform the classification on an image with an atmospheric correction (Dark Object subtraction), it returns an unclassified image.

The atmospheric correction rescaled the image on lower values, but I can't understand which scale factor I have to insert.


texts:
I have a dataset of several states and a table with many informations about them, they are already joined, what i want to do now is to create 5 classes about the average lifespan.

they have to be within the ranges of 37 - 50 for class 1, >50 - 60 for class 2, >60 - 70 for class 3, >70 - 75 for class 4 and >75 for class 5.

Now i have already tried the reclass by table tool but for some reason it wont allow me to use the needed fields. when i tried it via properties and then tried to type in 5 classes i couldnt do it as well.

im actually supposed to create a new field and do it by the field calculator in the table but since this field calculator is so different to the selection by attributes im also absolutely lost with that.

Does anybody know a way that might work for my problem?


texts: I need to apply the Ripley's K-function on a point pattern in QGIS 2.4 Chugiak, but I don't know how to do. I think I need to use some plugin or extension (maybe GRASS, R or something), but I haven't found neither clear nor structured solution about it.


texts: Creating map-based online questionnaire/survey?

texts:
I am planning to create map-based online questionnair/survey for my master dissertation on bicyle safety. What I want to do is, asking my potential respondants to mark their answers by adding pins/mark on Online map I create. For example, I ask them a question like, " where did you have a near-miss incident while you are on a bicycle?", then respondant can answer that question by locating pins/mark on where they had nearmiss.

My initial investigaton suggested that I can achieve this with Google map API, ArcGIS online or Crowdmap. I have seen some crowd-sourced maps with a lot of pins located by many people. But I am not sure if I can;

download all these marked locations in csv/text file with xy coordinate
make my survey map accessable only for my "group"

Can anybody suggest me what can I use (Google API, crowdmap, ArcGIS online or else) to achieve above 1 and 2 without any or with little programming?


texts: Transforming coordinates from one projection to another

texts:
var bounds = config.init_center.transform( new OpenLayers.Projection("EPSG:31466"), new OpenLayers.Projection("EPSG:3857"));

    console.log(config.init_center);
    console.log(bounds);

I have coordinates from a map that uses projection epsg:31466 and want to show then on a map with coordinates epsg:3857.

I would expect that the output of the two console logs given above would turn out to be different, since a transformation took place. However, they do NOT differ. They are exactly the same.

What is wrong with my assumption?


texts:
I have a dataset of pairs of parallel linear features (railway tracks) and I  need to simplify each pair into its central/median line. This could be achieved by computing a skeleton as suggested in this question but I am looking for a simpler solution: a line being the "average" of the two input lines would be a fair approximation. For example:

Would such transformation be implemented in any mainstream GIS software or library?


texts: Find extent of data inside many polygons

texts:
This is similar to Calculate polygon area inside extent

I have a few different polygon datasets which describe things like soil type, bushfire risk etc. For each property boundary we want to see how much of each soil type etc there is.

I tried Tabulate Interesection

but it gives me percentages over 100%  . The 3 in blue are the same property.

Any suggestions?


texts: How to increase performance for nearest neighbour query using a multipoint? (Boost Geometry)

texts:
The following function gets the nearest point in an rtree to a multipoint geometry. I'd like the performance of the function to be better. Is there any way to make it faster?

TimingLink* JourneyPattern::getClosestFromTimingLink(multipoint g,int maxFromSequence) {
    std::vector<LinkPair> result_n;
    linkFromTree.query(bgi::nearest(g, 1)  &&
        bgi::satisfies([maxFromSequence](LinkPair const& v) {
            return ((v.second)->fromSequence < maxFromSequence);
        }), std::back_inserter(result_n));
    if (!result_n.empty())
        return result_n[0].second;
    else
        return NULL;
}

The purpose of the function is to find a point in the rtree linkFromTree that is closest to the multipoint g and whose fromSequence is less than maxFromSequence. linkFromTree is an rtree created using the packing algorithm and it usually contains about 100 points. The multipoint g will also contain about 100 points.

The slow part seems to be using a multipoint (g) to do the nearest neighbour query, as I have a similar function that uses a point instead and it is much faster (average 20x faster).


texts: How to calculate the distance between a line and a point (using lat/lon coordinates)

texts:
In my Android app, I need to calculate the distance from a point (the device's location, which is a lat/lon coordinate) to a line (defined by two lat/lon coordinates), and get the result in meters.

I've taken a look at this page and the Java-code there, but I realize that I can't use that code (at least not with the lat/lon coordinates) because it will give me the distance/difference in degrees between the point and the line, and not in meters, and if I'm guessing right, it isn't possible to convert the difference in degrees to meters, right?

Therefore, I wonder how I can do this in another way, to get the distance in meters between a line and a point, both defined with lat/lon values.

EDIT: I took another aproach to the problem. Now I instead have a method which finds the point on the line which is closest to the device location.
This formula doesn't work exactly as intended when feeding it with lat/lon coorinates however. When the device location is on the line, it points right, but as the device movest away from the line, the error just gets bigger and bigger, until it stops at the end of the line.

Here is the function (based on the code in the first alternative of this answer):

/**
 * Get the closest point on a line.
 *
 * @param a The first point of the line.
 * @param b The second point of the line.
 * @param p The point to start at.
 *
 * @return The point at the line a->b which is closest to the point p.
 */
public static Location closestPointOnLine(Location a, Location b, Location p)
{
    // Store the vector a->p
    Vector a_to_p = new Vector(a, p);
    // Store the vector a->b
    Vector a_to_b = new Vector(a, b);

    // Find the square magnitude of a->b
    double squareMagnitude = square(a_to_b.x) + square(a_to_b.y);
    // Calculate the dot product of a->p and a->b
    double atp_dot_atb = a_to_p.x * a_to_b.x + a_to_p.y * a_to_b.y;

    // Calculate the normalized distance from a to the closest point
    double t = clamp(atp_dot_atb / squareMagnitude, 0d, 1d);

    // Create a new Location to store the result values in
    Location location = new Location(LocationManager.PASSIVE_PROVIDER);
    // Set the result coordinate values
    location.setLatitude(a.getLatitude() + a_to_b.x * t);
    location.setLongitude(a.getLongitude() + a_to_b.y * t);

    // Return the location
    return location;
}

The Vector class is just a class which holds two double values and represent a vector.
The clamp(double v, double min, double max) method clamps a value (v) between a minimum (min) and a maximum (max) value, in this case 0 and 1.
The square(double v) method calculates the square of a value (simply v * v).

I would be happy if someone could help me with either modifying this code a bit to return a correct point for WGS84/lat-lon coordinates, or suggest another way of doing it.
The result point should be as accurate as possible, lets say, less than 5 meters error.


texts: DB Manager returns NULL for GeoPackage layers

texts:
When I try to do some spatial query in the DB Manager (as Virtual Layer Provider) or in the Virtual Layer Window (Layer > Add new Layer > Create Virtual Layer...) it returns a NULL array, as it doesn't understand the GEOM field.
This happens with GPKG and KML. Although it doesn't happen with SHP. I tried with a KML layer (original), then I converted to GPKG and SHP. When queried in the BD Manager in GeoPackage Provider, it worked.
Is this an issue or exists an explanation?
I searched about it, but I didn't understand the conclusion.
DISCUSSION:
Geopackage - select on geometry column returns null ;
https://github.com/qgis/QGIS/issues/27914


texts: How can you join attributes from polygons to points which are within the polygon feature using PostgreSQL and PostGIS?

texts:
Possible Duplicate:
  How to perform a spatial join of point and polygon layers in PostGIS?

I don't know how to solve my actual problem.

I have polygons and points (some points are located in the polygons) and I want to add the polygon attributes to the points when they are inside the polygons.

Does anybody has an idea how to do this?


texts: How to coincide two geometries which has difference in positioning?

texts: I have two vector layers one representing countries and other representing states. I have states which is slightly dis positioned so that it is not exactly not coinciding with the country.Algorithm help is required .How can it be done in QGIS.Please find the image below for better explanation.The yellow should coincide with the outer linestring of green.


texts: What is the expression for "Current Geodatabase" to be used in ModelBuilder?

texts:
I have features in different geodatabases. When I process the features inside these geodatabases from Model Builder,and provide output names for them, I want them to be something like "CurrentGeodabaseName_%name%".

I want to know what is the expression for including Current Geodatabase Name in the output file.


texts:
I'm trying to install the rasterio python package using the conda package manager. But when I try to do this in a new clean conda environment, I get an "Found conflicts!" error with no information on what the conflicts are.

conda create --name newenv1 --channel conda-forge rasterio
Collecting package metadata (current_repodata.json): done
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.
Collecting package metadata (repodata.json): done
Solving environment: /
Found conflicts! Looking for incompatible packages.
This can take several minutes.  Press CTRL-C to abort.
                                                                                -failed

UnsatisfiableError:
Note that strict channel priority may have removed packages required for satisfiability.

I'm wondering whether this has to do with GDAL being previously installed on my computer, although as far as I can tell I have successfully removed it. Otherwise, I have no idea what is causing the problem. Does anyone know how to fix this issue?

Edit: I was able to successfully install rasterio from the pkgs/main channel, instead of the conda-forge channel that I had been using.

conda config --add channel pkgs/main
conda create --name newenv1 rasterio


texts: Merge attributes of points and lines using QGIS

texts:
I have a point vector-layer and a line vector layer. All points of the point layer are positioned on the lines (railway lines) of the line layer. I would like to add the attributes of the line layer to the appropriate data record of each point in the point layer by a geo-calculation (analogous to a geo-calculation of points within polygons).

How can this be accomplished?


texts: Problem with the layers display KML stored in GeoServer

texts:
I have a problem with the layers display KML stored in GeoServer. Google Earth shows me the study area, but the layer doesn't appear. However, I can see the layer when I "click" in OpenLayers.

Does anyone know what happens?


texts: GRASS shell commands not working in QGIS

texts:
Running GRASS 6.4.3 within QGIS 2.4.0 on Mac OS 10.10 Yosemite. I'm trying to run v.net.iso within the GRASS shell, but the following error is returned every time:

bash: v.net.iso: command not found

I did this successfully two weeks ago before I upgraded to Yosemite, but have not been able to get it work since. I've tried a number of other commands in the shell module, and I get the same error each time.


texts: Query Feature using OpenLayers + MapServer

texts:
I'm trying to query a feature in OpenLayers, but it's returning just an image highlighted that is not possible to make zoom, pan and drag.

My intention is to query a feature on the map like street name, city and country. Is there any function in OpenLayers?

I'm using this code:

var query = new OpenLayers.Layer.WMS( "Utbs",
    "http://localhost/cgi-bin/wms.exe",
    {
    srs: 'EPSG:4326',
    width: '800',
    styles: '',
    height: '600',
    layers: 'layer_name',
    mode: 'itemquerymap',
    qstring: "(value = '03')",
    qlayer: 'layer_name',
    mapxy: 'shape'
    },
    {
    singleTile: true
    }
);

map.addLayers([utbs_query]);
map.zoomToMaxExtent();


texts:
I am trying to convert coordinates (and polygons from com.vividsolutions.jts.geom) in Lambert72 (EPSG:31370) to WGS84 (or the other way around). I'm using the org.geotools.referencing CRS and MathTransform.
To check I used an existing point from an official government database.
As you can see, the transformed coordinates do not match but differ by up to 70 meters.

I'm using Java Oracle 8 and GeoTools 14 (with gt-epsg-wkt-14.0) in Eclipse.

Am I not using the MathTransform or the coordinate systems correctly?

Code:

public static Geometry convertWGSToLambert(Geometry g) {
    CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:4326");
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:31370");
    MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
    return JTS.transform(g, transform);
}

public static Geometry convertLambertToWSG(Geometry g){
    CoordinateReferenceSystem sourceCRS = CRS.decode("EPSG:31370");
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
    MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
    return JTS.transform(g, transform);
}

public static void main(String[] args) {
    //Point:
    //1. Lambert72:
    double Lx = 148349;
    double Ly = 170851;
    Geometry lambert = geometryFactory.createPoint(new Coordinate(Lx, Ly));
    // 2. WGS84:
    double wgsX = 50 + 50.0 / 60 + 53.0 / 3600;
    double wgsY = 4 + 20.0 / 60 + 43.0 / 3600;
    Geometry wgs84 = geometryFactory.createPoint(new Coordinate(wgsY, wgsX));

    // conversion tests
    System.out.println("Conversion test case");
    System.out.println("Lambert->WSG: " + convertLambertToWSG(lambert));
    System.out.println("    instead of      (" + wgsY + "\t" + wgsX + ")");
    System.out.println("wgs84->Lambert:" + convertWGSToLambert(wgs84));
    System.out.println("    instead of      (" + Lx + "\t" + Ly + ")");
    System.out.println("\n");
    System.out.println("\"Dubbele\" transformatie:");
    System.out.println(convertWGSToLambert(convertLambertToWSG(lambert)));
    System.out.println("instead of (" + Lx + "\t" + Ly + ")");
}

Output:

Conversion test case
Lambert->WSG: POINT (4.344261536168711 50.847806380256316)
    instead of      (4.345277777777778 50.84805555555556)
wgs84->Lambert:POINT (148420.57570480864 170878.70642028935)
    instead of       (148349.0           170851.0)

"Dubbele" transformatie:
POINT      (148348.99959955606 170851.00944651756)
instead of (148349.0           170851.0)

CRS WKT's:

GEOGCS["WGS 84",
  DATUM["World Geodetic System 1984",
    SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
  PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
  UNIT["degree", 0.017453292519943295],
  AXIS["Geodetic longitude", EAST],
  AXIS["Geodetic latitude", NORTH],
  AUTHORITY["EPSG","4326"]]

PROJCS["Belge 1972 / Belgian Lambert 72",
  GEOGCS["Belge 1972",
    DATUM["Reseau National Belge 1972",
      SPHEROID["International 1924", 6378388.0, 297.0, AUTHORITY["EPSG","7022"]],
      TOWGS84[-106.8686, 52.2978, -103.7239, 0.3366, 0.457, -1.8422, -1.2747],
      AUTHORITY["EPSG","6313"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH],
    AUTHORITY["EPSG","4313"]],
  PROJECTION["Lambert_Conformal_Conic_2SP"],
  PARAMETER["central_meridian", 4.367486666666666],
  PARAMETER["latitude_of_origin", 90.0],
  PARAMETER["standard_parallel_1", 51.166667233333335],
  PARAMETER["false_easting", 150000.013],
  PARAMETER["false_northing", 5400088.438],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["standard_parallel_2", 49.833333900000014],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","31370"]]


texts: Free good quality small scale (1:70,000) global coastline datasets?

texts:
I've been searching for a few days for the most accurate world coastline datasets.. Does anyone know where l can get coastlines at scales of 1:250 000, 1:2 000 000 and 1: 70 000. I've tried using the NOAA NGDC extractor but my 1:70 000 medium resolution coastlines only cover the US and the rest of the world is omitted.

I already have 1: 1 000 000(VMap0). Please point me in the right direction of where l can easily download good quality large-scale global coastline datasets.

Any help will be greatly appreciated


texts: In Qgis 2 I used to get high resolution geopdfs form the Georefexport plugin. In Qgis 3 i can get geopdfs by using the print composer/layout. However, I can not find a way to get high resolution geopdfs of many megabytes, so that when zooming into the geopdf, there is a lot of map detail, or "depth". How can I export high resolution geopdfs using the GUI in Qgis 3? Im on Qgis 3.4.2.


texts:
im trying to create a route displayer for ships. But im struggling with the edges of the map. When going from ie. 179* to 2* it displays a line across the entire map, instead of displaying it on the continuation of the map. (Map repeats)

I do see why it would do that. But i was hoping there was some support in openlayers for making this more or less automatic instead of breaking up the route.

Any ideas on how to fix this?

Expecation:

Current display:


texts: Reduce Processing Time for Update Cursor

texts:
I am updating a field for a feature class that is on an enterprise geodatabase in SQL server. This is an ArcSDE connection.

My goal is to search for null values, and convert these to 0. The script is currently run through python on ArcCatalog 10.3. This is the part of the script in question:

fc= r"C:\Users\migrate\AppData\Roaming\ESRI\Desktop10.3\ArcCatalog\OZ@DB_MAPS.sde\DB.OZ.BaseMap\DB.OZ.PlantCenter"
arcpy.MakeFeatureLayer_management (fc, "fc_temp")
with arcpy.da.UpdateCursor ("fc_temp", fieldID) as cursor:
    for row in cursor:
        if row[0] == None:
            row[0] = 0
            cursor.updateRow(row)

This dataset has about 20,000 entries, and the processing time is immense (over an hour at this point) for this part of the script. Is there a way to change the code to improve processing time for update cursors?


texts: Extract string between special characters - and / into a new field ArcGIS Field Calculator

texts:
Can I create a field calc expression in ArcGIS 10.3 to print a selection between special characters?

I would like to extract a selection of string from field [FolderPath] to a new field [Name]. The string is between a hyphen and a forward slash.

What I have in field [FolderPath]  example string:  264K - Name of Place/FeatureType

What I want in field [Name] target string: Name of Place


texts: Using ArcPy to produce lyr files pointing at geopackage?

texts:
I have hundreds of lyr files with a complex tree structure containing both layer groups and nested layers.
What I would like to achieve is to grab each lyr files, iterate through every single layer (also within a layer group), change the data source and the name of the dataset it points to and most importantly change the support from ArcSDE (underlying sql server DB) to a gpkg file that already contains the tables which the lyr will point to.

I produced a lyr file with all layers pointing to the wanted gpkg file and table manually, it works fine; I used it as a template for the other lyr files I want to change, so that I only have to change the source with the replaceDataSource method.

def update_from_template(self):
    # adds the template lyr to the temp mxd
    source_main_lyr = arcpy.mapping.Layer(self.template_lyr)
    arcpy.mapping.AddLayer(self.updt_df, source_main_lyr, "BOTTOM")
    for upd_lyr in self.layers_list:
        if not upd_lyr.isGroupLayer:
            # keeps track of the original datasetname
            datasetname = upd_lyr.datasetName
            # tracks down the source lyr (pointing to SDE)
            source_lyr = arcpy.mapping.ListLayers(self.updt_mxd, upd_lyr.name, self.updt_df)
            if len(source_lyr) > 0:
                # change the source lyr properties to the same of the template (pointing to GEOPACKAGE)
                arcpy.mapping.UpdateLayer(self.df, upd_lyr, source_lyr[0], False)
                # changes the name of the source back to its original
                newName = datasetname.replace("SDE_LIVE.DBO.GLOBAL_", "main.")
                # I use SDE again here (there is no GEOPACKAGE OPTION) I tried also NONE but I get the same
                # results
                servProp = upd_lyr.replaceDataSource(self.repointed_GPKG_location, "SDE_WORKSPACE", newName, False)
                upd_lyr.DataSource = servProp

But result of the lyr file are broken links, pointing to the old (template) dataset. I imagine what I am trying to do is not immediatly supported by ESRI but I was wondering there is anyone out there who had tried this, and maybe give advice. I am using ArcGIS desktop 10.4 and arcpy python 2.x but could update if required.


texts: Set image encoding default to PNG in QGIS WMS connections

texts:
Some WMS connections in QGIS offer a variety of image encodings and sometimes from am list [PNG, PNG8, JPEG, GIF, TIFF, SVG] JPEG ist selected per default, like shown below, sometimes PNG is selected per default.

Note that the default is selected any time one connects to the WMS an cannot be overridden by selecting another option.

Example #1: default JPEG

Example #2:default PNG

The default image encoding seems not defined in the WMS' capabilities, nor am I able to define it when adding a WMS connection manually or via PyQGIS like

wms = {'Digitales Orthophoto 80cm (LDBV)':[('', '', ''),(7, False, False, False, False, '', False, 'https://geoservices.bayern.de/wms/v2/ogc_dop80_oa.cgi')]}
set = QgsSettings()
set.beginGroup('qgis/connections-wms')
wms_exist = set.childGroups()
set.endGroup()
wms_cons = wms.keys()
for w in wms_cons:
    if not w in wms_exist:
        auth, params = wms[w]
        authid, user, pwd = auth
        dpi, iaxis, igfi, igmu, invaxis, ref, smooth, url = params
        wms_auth_path = 'qgis/WMS/{0}/'.format(w)
        set.setValue(wms_auth_path + 'authcfg', auth)
        set.setValue(wms_auth_path + 'password', pwd)
        set.setValue(wms_auth_path + 'username', user)
        wms_conf_path = 'qgis/connections-wms/{0}/'.format(w)
        set.setValue(wms_conf_path + 'dpiMode', dpi)
        set.setValue(wms_conf_path + 'ignoreAxisOrientation', iaxis)
        set.setValue(wms_conf_path + 'ignoreGetFeatureInfoURI', igfi)
        set.setValue(wms_conf_path + 'ignoreGetMapURI', igmu)
        set.setValue(wms_conf_path + 'invertAxisOrientation', invaxis)
        set.setValue(wms_conf_path + 'referer', ref)
        set.setValue(wms_conf_path + 'smoothPixmapTransform', smooth)
        set.setValue(wms_conf_path + 'url', url)

Since JPEG and PNG are available in both examples and the default values differ, the following question arises:

Q: How does QGIS determine the default image encoding and (if at all) how am I able to force QGIS to select PNG as default if available?


texts: Adding vector layer over WMS using OpenLayers?

texts:
I'm trying to add vector layer to my map from a GeoJSON file stored in the same directory as my maps index.html. However, it won't work.

Is there something wrong with my syntax?

My code:

<!doctype html>

    <link rel="stylesheet" href="http://openlayers.org/en/v3.0.0/css/ol.css" type="text/css">

    <title>Map of San Diego</title>

My Map


texts: NetCDF to Raster: Dealing with Memory Leak

texts:
I have a short script that I am using to process roughly 5,000 NetCDF files so that I can extract a single raster layer. This is not difficult, and I have a short script that does it as follows:

x_dimension = "lon"
y_dimension = "lat"
band_dimension = ""
dimension = "time"
valueSelectionMethod = "BY_VALUE"

outLoc = "C:/Out_location"
import_path = 'C:/import_location'

for root,dirs,files, in os.walk(os.path.abspath(import_path)):
    for f in files:
        outfile = os.path.join(root,f)
        write_file = ('%s' %outLoc)+f +".img"
        layer = 'name'
        arcpy.MakeNetCDFRasterLayer_md(outfile, variable, x_dimension, y_dimension, layer,'#', '#', valueSelectionMethod)
        arcpy.CopyRaster_management(layer,write_file, "#", "#", "#", "NONE", "NONE", "")
        del(write_file)
        arcpy.Delete_management(layer)

print 'Done'

This script works fine, but it has an aggregious memory leak that only allows it to process roughly 125 files at a time.

What is happening is that either function being used to convert my netCDF to raster automatically caches the raster it converts into memory and hence accumulates a raster with each loop that is run (which I doubt since I am using a temporary variable for that raster name that is overwritten with each loop); or that the function that creates a raster copy is caching the raster it creates for, for some reason, in memory.

Hence, I am trying to configure it such that my memory leak is solved so that I can process the entirety of my 5000 files at once.

I am beginner at using python, but from what I understand a way to halt a memory leak via using subprocesses, but quite frankly, even after my research, I am not sure how to do this.

Things I have tried:

I have tried deleting the temporary file that is created by raster creation tool but nothing happens
I have input a string that is taken as the name of the output raster for the NetCDF to raster tool that is input into the Raster copy tool and is overwritten every loop

I have tried deleting the output raster from the copy raster function in every loop which eliminated the leak, except it deleted it from my save location as well which indicates that this is a linked to raster creation function.

I am at a loss, hence any expedited solutions would be great!


texts: Is it possible to wrap text in column headers in print composer tables using QGIS 3.x?

texts:
Is it possible to wrap text in the header of an attribute table in print composer in QGIS 3.x?  There are options for wrapping text in the cells within the body, but they don't seem to work for the headers.
Other questions on stack exchange indicate it was not possible in previous versions, but all of the questions/answers are several years old, and there has been significant progress since then.  (This question, for example: Is it possible to wrap the heading of a column of an attribut table within the print composer?)
I'm using version 3.8.3 (Zanzibar)

I've tried both adding \n to the attribute heading, and setting the "wrap text on" box to * and adding that to the header.  Both of these approaches work within the body of the attribute tables.


texts: Keep one distinct overlapping feature in a shapefile

texts:
In the attached image you can see that there are features in the shapefile that have the same geographical extent, I would like to find a way to only keep one of them and remove the rest.

I would like to keep the one with the lowest value on the attribute Year. Using the remove duplicate geometries tool, I can't find any way to just keep the feature with the lowest value on the attribute Year.

Is this possible automate with any opensource tool like QGIS or GDAL? There are over 10 000 features so to do it manually is a bit too time consuming.


texts:
I was looking at the GDAL metadata of a USGS GeoTiff quad (o48118e8) with the GdalToTilesWin program. It dumped the following info:

Name o48118e8.tif
Size X/Y 5315/6780
ResolutionX/Y 2.4384/2.4384
Extent
North 5388204.0171
West 350547.3487
South 5371671.6651
East 363507.4447
It is not Geodesic WGS84(EPSG:4326)
PROJCS[
    "NAD27 / UTM zone 11N",
    GEOGCS[
        "NAD27",
        DATUM["North_American_Datum_1927",
        SPHEROID[
            "Clarke 1866",
            6378206.4,
            294.9786982139006,
            AUTHORITY["EPSG","7008"]
        ],
        AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4267"]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT[
        "metre",
        1,
        AUTHORITY["EPSG","9001"]
    ],
    AUTHORITY["EPSG","26711"]
]

Could someone explain how to understand extent? For instance, does North 5388204.0171 mean 5388204.0171 meters north of the NAD27?

Or maybe the Unit.Metre has something to do with it?

How do these 'extents' convert to lat/long (if at all?)

Also, out of curiosity, I presume the SizeX/Y is pixels, but what is the ResolutionX/Y?


texts:
EDIT: Link to repo containing the files. Open the index.html file in the imports and exports folders.

I can get everything except featurecClick to work.

My JS:

window.onload = function(){
    var cartoDbTableName = 'sipri_import_export_map_1950_2014';
    var domId = 'map';
    var mapStyle = document.getElementsByClassName('map-style');
    var lat = 0;
    var lon = 0;
    var zoomLvl = 2;
    var options = {
        center: [lat,lon],
        zoom: zoomLvl
    };
    var mapObject = new L.Map(domId,options);
    var layerSource = {
        user_name: 'chrismp',
        type: 'cartodb',
        sublayers: [
            {
                sql: "SELECT * FROM "+cartoDbTableName+" WHERE (gwsyear <= 1950 AND gwsyear > 0)",
                cartocss: mapStyle[0].innerHTML
            }
        ]
    };

    L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png')
        .addTo(mapObject);

    cartodb.createLayer(mapObject,layerSource)
        .addTo(mapObject)
        .on('done',function(layer){
                layer.getSubLayer(0).on('featureClick', function(e, latlng, pos, data, subLayerIndex) {
                  console.log(e, latlng, pos, data, subLayerIndex);
                }).on('error',function(err){
                    console.log('featureClick error: '+err);
                });
        }).on('error',function(err){
            console.log(err);
        });
};

The map loads, but I can't actually click any country. What gives?


texts:
Since a shapefile can only contain a single geometry type (e.g. point, polygon, etc), how could I persist both?

If I must create a shapefile for each geometry type, how can I pass them to a colleague to load?  In kml I can just pass the single kml to load.  Is there a shapefile manifest equivalent that I am missing?


texts: Can I export a pdf that renders my layers differently at different zoom levels?

texts: For example, say I have some streets I would like to show up, and streets are only set to appear at certain zoom level, well that layer does not get exported.  There must be a way to do this for pdf, am I right?


texts: Use Only First "x" Characters from Table in qgis:selectbyexpression

texts:
I have a list of street names in this format:

STRNAME
Elsinore
Neil Young
Penny
Fall

And a shapefile ("input_streets") with street names in this format:

STREET
Elsinore Av
Neil Young St
Penny Lane
Fall Dr

I want to select streets in the "input_streets" shapefile that are also found in the 1st list, but you can see my challenge, the first list doesn't have the street suffixes, the second list does.  This means a simple expression like this won't work:

expression = '"STREET" = '+ "'"+ strname"'"

...where "strname" is just the "STRNAME" values from the first list.  The suffixes in "STREET" values 'get in the way'.

I need a way to give QGIS an expression that can compare ONLY the Street Names (their Prefixes, Ignoring the Suffixes, whatever they happen to be: St, Lane, Rd, etc.)

Something like:

street = STREET[0:firstspace]

But can this be done?  Can one refer to the column heading of a table with RESTRICTIONS on HOW MUCH of the column value to USE in an expression? (ie. "give me the value in STREET but stop at the first space")

Example:

street = STREET[0:firstspace]
expression = 'street = '+ "'"+ strname"'"

Or can one only refer to the column heading in an expression, and whatever values happen to be in that column, you're stuck with?


texts: Finding the "square" bounding boxes of a country

texts:
I want to find the square bounding boxes that cover a country, of different sizes.

First off, I know that they won't really be square, but just going by the numbers they can be said to be.

The main challenge I have is I don't know where get the data of the points of the countries and/or the borders.  After that I think I can calculate the boxes by iterating over them.

Let's say I want do this for the continental US and with boxes of size 50 (lat, lon), then I can easily do this manually.  The first box (starting from NW) would be [[50, -150], [0, -100]], and the second is [[50, -150], [0, -100]], and I'm done.

Now let's say the size is 25, and now the first box is [[50, -125], [25, -100]].  I will end up with more boxes but a more defined area and overall smaller.  I guess another term could be "search grid" calculations.

Is there any publicly available datasets or APIs that I can use to determine if a bounding box or a point belongs to a country?


texts: How to resample polygon vertices so that they are evenly spaced in Arc 10

texts:
I have digitized the outline of many river channels (by both manually picking points and using streaming mode), but now I need the vertices to be evenly spaced (for input to another program). I tried to use the "densify" tool which was partially successful, but this adds additional vertices at the specified distance whereas I wanted to also simplify the sections that were streamed and have too many vertices.

Has anyone any suggestions for a tool I might try? I have no experience in writing scripts, but may be able to modify another. Many thanks in advance.


texts:
I have three fields of output vector layer in Qgis script:

##output_vector_first=output vector
##output_vector_second=output vector
##output_vector_third=output vector

When user fill first field I want recreate his path and use it for the next two fields.

For example:

User fill the first field: C:Users/myPath/first.shp
Then the next two fields will be autofilled with: C:Users/myPath/second.shp and C:Users/myPath/third.shp

Is it possible dynamically change field values in PyQgis script?


texts: Forcing a line layer's geometry to conform to another in QGIS

texts:
I have 2 line layers in QGIS representing rivers, type A and type B. I want to be able to merge them into a single rivers layer and populate an attribute telling me if they are of type A, type B or both. The resulting layer cannot have duplicate geometries.

The problem is that the layers were digitized using different basemaps and were then simplified using different thresholds. Therefore, the layers don't line up and one layer is more geometrically correct than the other.

I've included a diagram to illustrate what I mean. I need to keep the geometry of the red layer and snap the tributary shown in the purple box to it. The type attribute in the purple box should be "B", "A" in the green box and "Both" for the rest.

I've tried to merge the layers then use v.clean but the resulting geometry is somewhere between the 2 layers.


texts: Bringing Alias properties with Feature Class to Feature Class Tool?

texts: I am using arcpy.FeatureClassToFeatureClass_conversion(fc, outWorkspace, fc) to move many FC's to a new location.  I have found that it is only bringing the FC Name property and imposing that on the Alias property(see below).


texts: What parameters do azimuthal and transverse cylindrical projection in secant case have as proj4 string

texts:
conical equal-area projection, secant case has e.g.

+proj=aea +lat_1=... +lat_2=... +lon_0=... +lat_0=... +x_0=... +y_0=... etc...

I would ask if there is such thing, for transverse cylindrical equal-area, secant case?
e.g:

+proj=tcea +lon_1=... +lon_2=... +lon_0=... +lat_0=... +x_0=... +y_0=... etc...

and what is for the azimuthal equal area, secant case, I know just this for tangent case

+proj=laea +lat_0=... +lon_0=... +x_0=... +y_0=... etc...

what should one add to get secant case not only tangent at a point?

Edit:


texts:
When trying the fmask algorithm, QGIS crashes. This does not happen when using the blue band algorithm. I have used some versions but the result is the same.
I do not know if the problem is due to a missing python library.  I have tried to run the fmask algorithm from python but got no results.


texts: How are the GUID values for AddIn extensions created?

texts:
The AddInID GUID is straightforward, but when an extension is defined inside the AddIn,  it is also assigned a GUID value when it is enabled for the first time.

How (or where) is this GUID value obtained?

I've observed that it is the same across a couple of test machines.  Is this always the case?  If not, what would cause them to be different?

Example:

List of extensions before in the registry(HKEY_CURRENT_USER\Software\ESRI\Desktop10.x\ArcMap\Extensions) before enabling an extension on a freshly installed AddIn:

After the extension is enabled in the ArcMap Extensions dialog window:

The extension GUID 90814AD9-4737-553A-DF7F-E7039E4B3AF7 is now present in the extensions registry. This is not the same value as the AddIn ID.

I've searched the current user profile and the AddIn source code for this value with no results.  Where does it come from?


texts: gdal_merge not working in GDAL 201 windows install

texts:
I realise this looks like a duplicate of https://stackoverflow.com/questions/31174320/gdal-merge-py-not-working-after-gdal-and-it-python-bindings-have-been-installed, but I think the problem is a bit different.

I have installed GDAL201 MSVC 2013 win64, Python 3.4.4 and the corresponding GDAL-Python bindings on a windows 7 64 machine.
GDAL and Python both seem to be correctly installed. I have also set a system variable to the correct python location and I was able import GDAL and ogr into python.

However, when I run the command gdal_merge from the GDAL prompt I get the error: 'gdal_merge' is not recognzed as an internal or external command, operable program or batch file.
Gdal_merge.py gives me a "gdal_merge.py is not a valid Win32 application." popup box and an "access denied" error at the GDAL prompt.
Other functions such as gdalwarp seem to work fine since when I type gdalwarp at the prompt I get a usage message.

EDIT:
This is the original command from the script I was trying to use:

gdal_merge.py -o ..\merged.tif -of GTiff -init 255 -co "TILED=YES" -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" --optfile tiff_list.txt