texts:
I'm newbie and have question about string builder to show the arcgis layers.
(System: ArcGis 10.2, VisualStudio2012)
I don't succeed to show the loaded layers in a line of the message box.

code:

Imports ESRI.ArcGIS.Geodatabase
Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.ArcMapUI
Imports System.Text

Public Class LayerLoop2_Button1
  Inherits ESRI.ArcGIS.Desktop.AddIns.Button

  Public Sub New()
  End Sub

  Protected Overrides Sub OnClick()
        AccessLayerViaMxDocument(My.ArcMap.Document)
        My.ArcMap.Application.CurrentTool = Nothing
  End Sub

    Shared Sub AccessLayerViaMxDocument(ByVal pMxDocument As IMxDocument)
        Dim pMap As IMap
        Dim pMaps As IMaps = pMxDocument.Maps
        Dim i As Integer

        For i = 0 To pMaps.Count - 1 Step i + 1
            pMap = pMaps.Item(i)
            Debug.WriteLine(pMap.Name)
            MsgBox(pMap.Name)

            Dim pEnumLayer As IEnumLayer = pMap.Layers(Nothing, True)
            pEnumLayer.Reset()
            Dim pLayer As ILayer = pEnumLayer.Next()
            Dim sb As System.Text.StringBuilder = New System.Text.StringBuilder()

            While Not pLayer Is Nothing

                Debug.WriteLine(pLayer.Name)
                ' Show MessageBox with the name of the layer
                pLayer = pEnumLayer.Next()

                sb.Append(pLayer.Name)
                MsgBox(sb.ToString())
            End While
        Next
    End Sub

Can you explain me please what you mean with console output? As output I should just have a message box. In the code above I tried already to use the string builder - but without success....
The first given code was this:

Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.ArcMapUI

Public Class AccessLayer

Inherits ESRI.ArcGIS.Desktop.AddIns.Button

Public Sub New()
End Sub

Protected Overrides Sub OnClick()

AccessLayerViaMxDocument(My.ArcMap.Document)
My.ArcMap.Application.CurrentTool = Nothing
End Sub

Shared Sub AccessLayerViaMxDocument(ByVal pMxDocument As IMxDocument)

Dim pMap As IMap
Dim pMaps As IMaps = pMxDocument.Maps

Dim i As Integer

For i = 0 To pMaps.Count - 1 Step i + 1
pMap = pMaps.Item(i)

Debug.WriteLine(pMap.Name)

MsgBox(pMap.Name)

Dim pEnumLayer As IEnumLayer = pMap.Layers(Nothing, True)
pEnumLayer.Reset()

Dim pLayer As ILayer = pEnumLayer.Next()
While Not pLayer Is Nothing
Debug.WriteLine(pLayer.Name)

MsgBox(pLayer.Name)

pLayer = pEnumLayer.Next()
End While

Next
End Sub

Protected Overrides Sub OnUpdate()
Enabled = My.ArcMap.Application IsNot Nothing
End Sub

End Class


texts: Deleting selected features using PyQGIS?

texts:
I have a shapefile with only one field DN. I would like to delete all features which does not have specific value of the field. Let's say the DN can be 1,2 or 3.

So far I was able to find out how to select Features and how to delete feature, but I did not succeed to combine them. Let's say I want only features with DN = 3 to remain:

Selection of the Features found here, however, I do not know if it does what it is supposed to do since I did not find the way how to check what is inside.

select = layer.getFeatures(QgsFeatureRequest().setFilterExpression(u'"DN"!=3'))

delete features found here

res = Ilayer.dataProvider().deleteFeatures()

As far as I know, the deleteFeatures() requires IDs of the features, but I have no idea how to get them from from select which is QgsFeatureIterator.


texts: How to create a symbol with arrows and an arc in QGIS

texts:
I'm using QGIS 2.18 and wondering if there's a fairly simple way to create a symbol with an arrow on the end of a ray, and an arc in between.

The thing that makes it complex is that the angle is not consistent.  Also, I always want the arc to travel from the first line drawn to the second line drawn so I can show angles greater than 180.  In my drawing point 1 is the left most arrow, point 2 the center, and point three the other arrow.  So the arc would travel from line 1-2 to 2-3.


texts:
I installed i.ann.maskrcnn in GRASS GIS 7.8.3
but, I got error message in the picture.

How can I download metadata file about this tool?


texts: Is there a way to use the project layout if using a plug in that creates batch maps through QGIS?

texts: There is a plugin for biological records which creates batch maps for each taxon but it bypasses the premade project layout with a legend, picture of a logo and the copyright for each of them. Is there a code to run the batch creation through the project layout?


texts:
Is there a reason that every script I have written for GRASS in python runs twice?  I start them using the GUI script launcher tool, and they all begin with:

#!/usr/bin/python
import grass.script as grass

So nothing fancy.  Still, every script always runs through twice, first one run through and then another, with the first run through locking the GUI until it finishes, and the second unlocking the GUI & displaying status update bars, etc.  I've confirmed that they're running twice both by using log files & with ps -e.  This happens both in GRASS 7.0 and in GRASS 6.4.2 (Ubuntu 12.04).  Am I missing something obvious?  If not I'll go ahead & submit a bug ticket.


texts:
There is a nice resource on unit testing for QGIS. There is also a post on this, Writing automated tests for QGIS plugins, which requires the use of docker and relies on QGIS 2 plugin - QGIS Tester.
The progress on QGIS 3 is related to two plugins: 1) Tester and 2) QGIS Core Tests.
How this all ties up with the final screen of QGIS Plugin Builder once we activate the Unit Tests checkbox? Where to start when we want to unit test the developed plugins?


texts: Using "IN" CQL function in Geoserver Cross-Layer Filtering

texts:
I've two layers "cities" (has the fields: fid, code, and geom) and "cities_details" (has the fields: id, city_code, and area).

I need to filter "cities" layer with a condition applied on the "cities_details" layer. I believe that Geoserver's Cross-layer Filtering extension can do that.

The following CQL filter runs successfully:

code = querySingle('cities_details', 'city_code','area = 5000')
It returns the city with area = 5000.

But the following CQL filter failed:

code IN (queryCollection('cities_details', 'city_code','area >= 5000'))

I tried also the following, but didn't work:

code IN (collectGeometries(queryCollection('cities_details', 'city_code','area >= 5000')))

In both cases the log tells:
org.geoserver.platform.ServiceException: Could not parse CQL filter list

Can anybody help me, how to use "IN" with the Cross-Layer Filtering


texts: Transparent background in a multiband raster using GeoServer

texts:
I have a raster with a black collar (background). I want to set those black pixels to transparent. The problem is that the values of the pixels have three bands (Red, Green and Blue) and I don't know how to define a rule in the .SLD to set those pixels which value is (0, 0, 0) to transparent.

Any idea on how to do this?

My question has been associated with another question. The answer to the other question resolves the problem for a single-band raster but in my case I have a multi-band raster, so the answer I need is different.


texts: Generate GCP table from Geo-referenced raster in QGIS

texts: I've Geo-referenced raster, Is there is any way to extract it's GCP table in QGIS ?


texts: Unable to permanently turn off time stamp display in QGIS Time Manager

texts:
Using QGIS 2.10.1 Time Manager, I wish to permanently turn off the time stamp display for one particular project.  In Time Manager Plugin > Settings I uncheck "Display frame start time on map".  The time display correctly disappears.  However, when I save the project file and then reopen it later, "Display frame start time on map" is now checked on, and the time stamp displays.

How do I permanently turn off the time display for this project file?


texts:
I want to do a conversion from WCG to UTM for my 172 Raster. I can't use normal batch because I can't re-arrange again the order of raster in the Table of Content.

Because the result is a random name in the same folder, so I need to do Raster Projection Conversion without changing the name and extension, only the output folder.

How to Batch Raster projection conversion without changing a name and extension only changing the output folder?

Some of Python code:

import arcpy, sys

InFolder = r"c:\start"
OutFolder = r"c:\end"
OutSR = r"WGS_1984_UTM_Zone_50S"

arcpy.env.workspace = InFolder

for Ras in arcpy.ListRasters():
    arcpy.AddMessage("Projecting " + Ras)
    arcpy.ProjectRaster_management (InFolder + "\\" + Ras, OutFolder + "\\" + Ras,OutSR)
arcpy.AddMessage("Projecting complete")


texts: GRASS and SAGA CentOS Linux release 7.1.1503 (Core)

texts:
I have a problem. I am developing a standalone app using QGIS but my script imports GRASS nad SAGA plugins and I have error to create the maps.

I have CentOS Linux release 7.1.1503 (Core). SAGA is not recognized.

Can you explain what I have to do?


texts: Finding vectors with wrong direction in river network

texts:
I have a shapefile containing network of rivers for the whole country. The description of the layer is saying that orientation of the vector should be the same as the river flow. This is true for most of the lines, but not for all of them.

I would like to find the vector with the wrong orientation and change it. I have already found a way with python how to change direction.

The problem is selecting all the broken elements: dataset is big(more then 300 000 elements), so that it was I cannot do it manually and looking for automatic solution.

I am working with qgis 2.14 or 2.15 and I am familiar with python.


texts: Exporting attribute data from vector data by selecting object in QGIS?

texts:
I'm using QGIS (2.4.0), and i need to export the attribute data of the object i select (shapefile) to another file ( Excel, txt,...it is yet to be determined depending on the solution).
I know how it is done manually, but this action has to be done over a thousand times, and it can only be done one object at a time...
I'm looking for a automated proces where just by selecting a object the attribute data is written on the fly into the other file.

Does anyone have a solution, or an easy workaround, for this?


texts:
I have a vector shape with many different sized polygons each with an other id and values. I want to assign a new table with Information  of each polygon containing the name of the majority (higher proportion) of adjacent polygons attribute name. Ergo some kind of polygon neighborhood analysis.

The query should look to the attribute of each bordering polygon and meassure the length of neighborhood polygon and id.

Does anyone have an idea or a tool how to accomplish this task?

I intend to use QGIS, or postgis. I know there is a tool in arcgis called neighborhood statistic.


texts: Saving project with data source path as relative in QGIS?

texts:
Is there an option in QGIS to save the project with the data source being a relative path?

I don't see the option here. QGIS Wroclaw1.7.3


texts: Using Asterisks in Python parser of ArcGIS Label Expression Builder?

texts:
I am trying to create a label expression that labels features that contain hike in their name (versus those that don't).

I have tried something like this:

def FindLabel ([Name]):
    if [Name] is "*hike*":
        return [Name]

But that eliminates all of my labels, instead of just labeling the hikes with hike in their name.

How would I do this in the label expression builder?


texts: Does Avenue equivalent of pointPosition exist in ArcGIS 10.2?

texts:
There were a couple of awesome methods in Avenue for ArcView:

point=polyline.Along(percentage)
percentage=polyLine.pointPosition(point)

Don't laugh, but one of the reasons I updated to 10 from 9, was new method for Geometry

point=polyline.positionAlongLine (value,{use_percentage})

that is equivalent of 1st Avenue method (wrong name in my opinion anyway). Does equivalent of the 2nd Avenue method exist in 10.2?


texts:
In ArcGIS 9.3.1 (and prior) whenever the user wants to add data to Arcmap from a specific folder (and not from the general drives, C:\ for example), He has to press "connect to folder" and map the specific folder. Other than being a nuisance, it doesn't allow the user to be flexible - if he suddenly needs data from a different folder he can't simply navigate to it, he has to map it first.

Does anybody know if it the same or was it fixed in ArcGIS 10?


texts:
I am a new GIS user, I have some problem can not be solved by myself. I get a shape file generated from an Excel file. at first, my projection is WGS_1984_UTM_Zone_48N some data are inside my country but some are out of the country. then I had selected and exported only the ones outside the country and reproject as WGS_1984_UTM_Zone_47N it is perfect. so my question is how could i display the two layers in the same projection( WGS_1984_UTM_Zone_48N ) perfectly. it would kind if anyone can help me.

I am using arc gis 9.3


texts: How get map extents for drawn polygon when final click is done?

texts:
I am using OpenLayers JS to draw a polygon and now I want map extents of the drawn polygon one final click (i.e. user finish drawing polygon or double click is done). For drawing the polygon I am using following code:

action = Ext.create('GeoExt.Action', {
        text: "Draw Poly",
        icon: 'Images/poly.png',
        control: new OpenLayers.Control.DrawFeature(
            vector, OpenLayers.Handler.Polygon
        ),
        map: map,
        // button options
        toggleGroup: "draw",
        allowDepress: false,
        tooltip: "Draw Polygon",
        // check item options
        group: "draw"
    });

This code is written using OpenLayers JS, Ext JS and GeoExt JS.

Can anyone please tell me how to get map extents for the drawn polygon?


texts: ModelBuilder summary statistics for one shapefile produced from clip iteration?

texts: I'm very new to ModelBuilder and trying to build a tool which clips 4 national data sets to a study area defined as a local authority. There is a points dataset containing population data, local authority boundaries and flood extents.  Ultimately the output I am after is a summary statistics showing the population within the flood extents for that study area.  I've managed to get the first part working where the user selects a local authority from a drop down list and the data is clipped to this study area through a feature class iterator.  This is a sub-model within my main model.  I then want to perform separate geoprocessing on each of the clipped feature classes, and this is the part I can't get working.  I need to work with the results of my sub-model, and produce summary statistics where the population points are within the flood extents.  However I can't seem to select the feature class to work with from the output values of my sub-model.  I've attached an image of my model for reference.


texts: Calculating band depth image for M3 data

texts: I am working with chandrayaan M3 hyper spectral data. I need to create a band depth image for identifying mineral abundances. How do i calculate the continuum reflectance value for producing a band depth image? Is that a constant value for the whole band or should i calculate the continuum reflectance for each pixel of the band. The formula given for calculating band depth is BD = 1- Reflectance at a specific wavelength(R1450)/Continuum reflectance at that wavelength(R1450i).


texts:
The following picture shows the location of the files I need for the Closest Facility Analysis.

What's wrong with my code?

import arcpy

arcpy.env.workspace = r"C:\Users\TL\Desktop\test\2.gdb"
inNetworkDataset = "C:\Users\TL\Desktop\test\2.gdb\network\network_ND"
inIncidents = r'G_Fac'
inFacilities = r'HighwayJunctions'
outGeodatabase = r"C:\Users\TL\Desktop\test\2.gdb"
outRoutes = "Routes"
outDirections = "Directions"
outClosestFacilities = "ClosestFacilities"
measurement_units = "Feet"

# Run FindClosestFacilities. Choose to find only the closest facility.
arcpy.na.FindClosestFacilities(inIncidents, inFacilities, measurement_units,
                                    inNetworkDataset, outGeodatabase, outRoutes,
                                    outDirections, outClosestFacilities,
                                    Number_of_Facilities_to_Find=1)

print "Script completed successfully"

The error I got:

Traceback (most recent call last):
  File "<editor selection>", line 14, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\na.py", line 2630, in FindClosestFacilities
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Network Dataset: Dataset C:\Users\TL\Desktop  est.gdb
etwork
etwork_ND does not exist or is not supported
Failed to execute (FindClosestFacilities).


texts: Ordering polygon centroid points in a row(s)

texts:
I have a image from a field, where I converted the plants to polygons, and then I found the centroid of those polygons (attached image). The numbers displaying on the image are the FID numbers.

I am looking for a way to either create new field on my attribute table showing the ID on the same order showing on the image below or to create a new ID field for each point (starting on the east side and moving to the west, for example).

I tried to sort the data set (around 35,000 points) by X and Y, but I does not work for me, and I assuming that is because I the points are just few inches apart.

My final goal is to measure the distance between each point (not the near distance), such as between 59-60, 60-50, 50-51, 51-33, and so on. Looking around on posts posted on this website, it seems that easier solution for my distance question is calculating them in Excel.

In order to do that I need to have the points in order, so I can calculate the distance to the next point.

I am very familiar with ArcGIS but I am not familiar with Python coding, but I can give a try if I need to.


texts:
I have a number of points with a "time" property. I need a procedure that, given a total time, results in a number of clusters in which the sum of the time property do not exceed the given total time.

Example below shows a clustering when the given total time is 15:
Blue cluster is 15
Green cluster is 15
Yellow cluster is 15
Cyano cluster is 14


texts: Publishing map created by qgis2web on CMS

texts:
I have a problem with the publication of a qgis2web export on a CMS (another than Wordpress). I have an error when I put .js/.ttf/.woff2 files:

these files are not accepted by the CMS.

Someone has already had this problem?


texts: FME doesn´t accept my ArcGIS 10.3.1 Installation as being current

texts:
When creating a lasd file with fme I get the following error although I have ArcGIS Desktop 10.3.1 installed

  FME has checked out an Esri license. The product checked out is
  'ArcInfo'

  LAS writer: ArcGIS 10.1 or newer is required to create an ArcGIS LAS
  Dataset.  No ArcGIS LAS Dataset will be created.  LAS files produced
  by the FME LAS writer have not been affected

Is there anyone who has defeated this error?


texts: Adding sequential codes to GPS waypoints

texts:
When using a GPS (a Garmin GPSmap 62 -- but our team has other similar equipment), when I want to add a new waypoint, the system displays a default number for it (e.g. 202). If I add the 202, the next point will be 203 and so on.

But, since the points I want to add belong to a project, I want to add a different code, intrinsic to the project (like ABC123, from project ABC, number 123). So I want to add the point ABC123, and configure the GPS in a way that the next suggested number will be ABC124, not 202 again (which causes an awful waste of time, after some points).

Does anybody know if this configuration exists in my device or perhaps in a different, more advanced model?


texts: Sensor Observation Service Filter Request

texts:
I am a student and currently developing a web application using Sensor Observation Service. In my Web-Client, I am sending request to SOS Server per sensor and getting back the observations.
The main functionality of Web-Client is to perform Spatial Selection on map and see the Observations of the selected region.

In this scenario, the region selected gives street/sensor Id's(which is fetched using WFS), and the client sends request to SOS for each street/sensor Id.
So, it means if I select 1000 street Id's, I have to send 1000 request to SOS, which will obviously slow down my application.

My question:

Is there any way, I can send filter request to SOS ? Like, I send all my sensor Id's to SOS and get the observations of all at once?
Or any other suggestion to improve my application behavior a part from SOS Filter Request?


texts: Merge and dissolve overlapping polygons?

texts:
I have a situation in which I have both network buffers and circular buffers around the same points (as separate feature classes), and I want to merge the buffers and dissolve into one polygon for each point.

The complication arises in that some of the points are close together so their buffers overlap; thus when I dissolve, both the network and circular buffers dissolve, but also the boundaries between buffers that are close together. But I need to maintain one row for each point feature.

Image shows the overlap between network and circular buffers, but also points with overlapping buffers.


texts: Easiest way to create historical street maps using QGIS

texts:
I want to make some historical street maps using QGIS. I understand how to upload a historical OS map from DigiMaps to QGIS, but I want to actually 'draw' the streets myself in order to have some detailed maps. The area I am dealing with was subject to a slum clearance in the 1960s and I am looking to map it how it was in the late-nineteenth and early-twentieth centuries.

Could anyone offer advice?

Thanks in advance!


texts: Converting file geodatabase with coded domains to PostgreSQL / PostGIS?

texts:
We have several file gdb's and would like to convert these to Postgres / PostGIS. We've used ogr2ogr to transfer these files; however, some of the fields are stored as coded domains and the fields are not decoded in the transfer. Coded domains are stored in a separate table.

For example, a column design quality could belong to the general domain quality. Quality could then be coded as '1=low, 2=medium, 3=high'. This is to enforce that only certain values are used, for instance in a survey. Ogr2ogr reads the values 1, 2, and 3, but not the accompanying text.

Would anybody know how to include these domains in the transfer from gdb to Postgres / PostGIS?

We could convert to file gdb's to an ArcSDE database if this is easier; however shapefiles are not preferred due to column truncation.

This question was posted earlier here on StackOverflow.


texts:
How to change attributes value when i move(change position) point?
I don't found some information about this.
Sorry for my english.


texts: How can I add a delete button to feature forms?

texts:
I have created some formulars in QTDesigner for capturing and manipulating data in QGIS. The data is stored in a postgis-database.

It would be helpful to have a delete-button in this formulars to delete the open feature. At the moment my formulars only contain the 'OK' and 'Cancel' Button.

How can I add a 'Delete'-Button?

Here my code. It's based on the example in the refered post. The delete() function is the missing part. Can I use something like currentItem,..?

from PyQt4.QtCore import *
from PyQt4.QtGui import *

id_flaecheField = None
id_fotoField = None
myDialog = None

def formOpen(dialog,layerid,featureid):
   global myDialog
   myDialog = dialog
   id_fotoField = dialog.findChild(QComboBox,"id_foto")
   id_flaecheField = dialog.findChild(QComboBox,"id_flaeche")
   buttonDelete= dialog.findChild(QPushButton,"pushButton")
   buttonDelete.clicked.connect(delete)

def delete():
# missing part of the delete function


texts: Create subsetted multiband image from multiple files using gdal or python?

texts: Given a time series of an area (single band), how you subset them and combine the subsets into a single multiband file (eg. vrt) using gdal or gdals' bindings for python?


texts: Grass 7 no longer works via the Processing Toolbox on Ubuntu

texts:
When I try to execute a Grass 7 algorithm via the Processing Toolbox in QGIS 2.10, 2.8.2 or 2.8.1, I get the following error:

  Missing dependency. This algorithm cannot be run :-(

QGIS erroneously claims Grass 7 is not installed, but it was working fine a few months ago. I'm on Ubuntu 14.04.

I installed Grass 7.0.0-1~exp2~trusty3 on my desktop and laptop using the UbuntuGIS Unstable PPA. I installed QGIS 2.10.0 (desktop) and 2.8.2 (laptop) using the software repositories published on the QGIS website.

In an effort to track down the problem, I downgraded my laptop from QGIS 2.8.2 to 2.8.1 using the UbuntuGIS Unstable PPA; the same PPA that provides Grass 7. However, I still get the same error. So I suspect problem is with Grass 7.0.0-1~exp2~trusty3 that was compiled at the end of March.

For a time a few months ago I had Grass 6 and 7 installed along side each other, I was able to access both via Processing in QGIS. Grass 6 is no longer available in the UbuntuGIS Unstable PPA, so downgrading Grass 7 may be impossible.

Does anyone have a fix for this or are you experiencing the same problem?


texts:
I'm trying to use the Vector > Analysis Tools > Sum Line Lengths in QGIS version 3.2.3-Bonn.  I'm inputting a ploygon (Multi Polygon) and a Lines (MultiLineString) shape files.  I'm hoping to get the count and length of lines in each square of the polygon grid.  The error I'm getting is;

  Algorithm 'Sum line lengths' starting… Input parameters: {
  'COUNT_FIELD' : 'COUNT', 'LEN_FIELD' : 'LENGTH', 'LINES' : 'X:/SAERI
  Contracts/CFL/Toothfish ms/Map/layers/LL_lines.shp', 'OUTPUT' :
  'memory:', 'POLYGONS' :
  'MultiPolygon?crs=EPSG:4326&field=fid:double(20)&field=left:double(23,15)&field=top:double(23,15)&field=right:double(23,15)&field=bottom:double(23,15)&field=id:long(10)&uid={411cf258-ceb3-4878-810c-8c14545a6a79}'
  }

  Traceback (most recent call last): File
  "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\algs\qgis\SumLines.py",
  line 119, in processAlgorithm
  QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(poly_source.sourceCrs(),
  context.transformContext())), feedback) Exception: unknown

  Execution failed after 0.02 seconds

Both layers are in the same projection (WGS84) and I've also tried UTM21S but get the same error.

Or is this a bug?

Or something wrong with my installation?


texts: Performance Monitor of ArcGIS for Desktop?

texts:
Is there a way/tool that could basically report what ArcMap 10.2 is doing?

I can't really explain it, so I will show you what I mean. Now mind you it doesn't have to be exactly like this.

ARCMAP Report
Query Basemap (150ms)
Draw Basemap (300ms)
Query Layer1 (100ms)
Draw Layer1 (200ms)
Query Layer2 (50ms)
Draw Layer2 (3000ms)

The point of this would show me that Layer2 is taking 3 seconds to draw, so maybe I could change the symbolism to make it faster.


texts: How to efficiently access files with GDAL from an S3 bucket using VSIS3?

texts:
So, GDAL has recently added a new feature that allows random reading of S3 bucket files. I am looking to crop GDAL images from multiple tiles of an image without having to download the whole file. I've only seen very sparse documentation on how to configure and access an S3 bucket though GDAL and am a little confused on how to begin? Would someone be kind enough to provide an extremely short example/tutorial on how one would go about setting the virtual filesystem for GDAL in order to accomplish this goal? Bonus pts if your solution allows it to be scripted via Python!

To clarify: We already have it done in Python. The issue with Python is that you have to download the whole image to operate it with it. The newest version of GDAL has support for mounting the S3 bucket so that if we need to say a crop a small portion of the image, we can operate directly on that smaller portion. Alas, as the feature only was released on the stable branch in January, I haven't found any documentation on it. So the solution should use the VSI3 system in the newest release of GDAL or otherwise smartly uses the system to prevent the user from needing to download the entire image to an EBS drive to operate on it.

That is to say the bounty will be awarded to answer that uses the VSI APIs found in the newest versions of GDAL so that the whole file does not need to be read into memory or disk. Also, we the buckets we use are not always public so many of the HTTP tricks being posted won't work in many of our situations.


texts: Manually changing/suppressing z-values in Spline with Barriers interpolation that are known to be erroneous?

texts:
I've included an image of a project I'm working on. There are two layers representing sub-surface geologic formations, and these were generated using elevation borehole data (so, point data), to interpolate each surface using Spline with Barriers. They will eventually be apart of a groundwater model.

The circled orange spots are problem areas where the interpolation appears to be incorrect. Is there a way to manually change these values myself so we can insert, say, an average instead?
(We are fairly certain these two layers do not "swap positions" in these spots.)

I haven't found much in the way of manually manipulating floating-point rasters after having generated them through interpolation, but I'm guessing there's something in the ArcToolbox that I'm just missing. I'm rather new to ArcMap.

I'm using ArcMap 10.3 Desktop with 3D Analyst and Spatial Analyst licenses.


texts: Find longitude intersection of given latitude at z=0 for a globe

texts:
I have a globe image that is centered at lat,lon. At this lat,lon, by definition x=0, y=0, z=1.

A line of latitude drawn on this map can either be completely visible (in which case it would be drawn as an oval like the 60°S line in the image below), or it can be completely hidden like the 60°N latitude in the picture (all points on the line are z < 0).

The third possibility is that a line of constant latitude has a start point at some longitude and an end point at some other longitude where z=0. How can I find these two longitudes?

For example if the center lat,lon is -50, -130 (between New Zealand and S. America), the line of latitude near 40°S starts near a longitude near Western Australia and ends in The Atlantic between S. America and Africa. But at what longitude does this line get to z=0 on both ends?

So in the picture below, the two red dots are on either side of the 40°S line. What are the two points of longitude? They will be at z=0, lat=-40 (in degrees) and long=???

In the picture below, gridlines are shown every 10 degrees with the center of the projection at Lat/Lon (-30,0). In this view, all latitude lines above 60N are invisible because they fall on the backside of the earth where z<0. All latitude lines below 60S are completely visible because the whole line in in an area where z>0.

The area between 60S and 60N has some of each latitude line visible. They end at the following longitudes (calculated by hand with trial and error and interpolation using excel):

60N   0.000
50N  46.523
40N  61.023
30N  70.529
20N  77.869
10N  84.157
Equ  90.000
10S  95.843
20S 102.130
30S 109.471
40S 118.977
50S 133.477
60S 180.000

I need a formula to calculate these values.

Given: A latitude for the center of the projection (-30 in this case).

Given: A line of latitude (20 for example)

Find: The number of degrees E/W that the line remains visible (z>=0): 77.869 in this example.


texts: QGIS equivalent of ArcMap's Selection by Location Have their Centroid In option

texts: I am currently using QGIS 2.18.15 and I have two sets of polygon layer data. I would like to use the Select by Location research tool in order to isolate the polygons that have their centers within the other separate layer. ArcMap has a specific "have their centroid in" option when selecting by location. I don't see that QGIS has this. Is there an equivalent?


texts: Markov chain to predict LULC in GEE - is that possible?

texts:
I was wondering if is possible to build a Markov chain model do predict future land cover with a collection of classified images.
Any tips?


texts:
I'd like to do a GeoJSON dump of a PostGIS table, but I want to export features (the geometry and the properties) not just the geometry. I've been digging into the ST_AsGeoJSON function but it looks like I can only get the coordinate information (which makes sense since it's a geometry function).

For example:

select st_asgeojson(the_geom) from street_centerline limit 1;

Returns:

{"type":"MultiLineString","coordinates":[[[-65.591776562805038,41.682190576167052],[-65.591998971284028,41.682082119060382],[-65.592001213509064,41.682081025737766],[-65.593689871787177,41.681257533373952],[-65.595415661879244,41.680415888937219],[-65.595440519465640,41.680403765889309],[-65.595603134242481,41.680324459445771]]]}

I'm curious if anyone knows of a simple way to get the property information as well. Has anyone written a pgsql2geojson script yet?


texts: Customising QGIS Map Composer default output directory and file name suffix?

texts: I am using QGIS for map making , I am adding multiple vector and raster layers and then finally for map making I am moving to the map composer. When I am exporting the map from Map composer , I want to export it to a particular directory and I want to add a suffix to the filename that user inputs in the Export dialog.


texts: Is Hot Spot analysis appropriate if you see hot spots in areas that you would expect to see them in?

texts: I am working with vessel position data on ArcGIS 10.0. My customer has requested density analysis maps. I was curious if hot spot analysis would be a more appropriate analysis? The study area is along a coastline of a country, and of course, there are going to be hot spots in choke points (port areas). Would hot spot analysis really show us any statistically significant areas that we wouldn't normally see? I understand that this completely depends on the dispersal of the data and it certainly would if you see an area that you wouldn't think to be a hot spot but was. I am having a hard time ensuring I am on the right track here instead of just providing a density map which everyone seems to want but that is entirely subjective.


texts: How to hold/collect values in an iterated modelbuilder

texts:
I'm using modelbuilder to do upstream traces from several points/flags in a feature class.  I'm iterating each point in the fc using the feature selection iterator so that I can calculate/transfer the inline variable/point ID to the trace results.  My question is, I have the output of the trace tool connected to the collect values tool which then outputs to the merge tool.  I want the merge tool to run only after all points have been iterated/traced and output into the collect values tool.  I thought this was the point of the collect values tool, that it would hold/collect all iterations and after the iterations were complete, only then would it continue to the next tool.  Currently the model iterates through all tools including the merge tool which then puts the error:

Failed to execute. Parameters are not valid.
ERROR 000400: Duplicate inputs are not allowed
WARNING 000725: Output Dataset: Dataset E:\Workspace\ScratchWorkspace.gdb\Merged_SSO_Assets already exists.
Failed to execute (Merge).

I'm guessing this error is due to a merge being performed at each iteration because when I run each tool individually in sequence, the model runs correctly and w/out errors.  I've tried using preconditions but they have not resolved the issue, I could be using the wrong preconditions.


texts:
I have a simple SLD style using in GeoServer v2.3.0 and testing in uDig.  I have a interstate symbol that I do not want the label symbol to follow the line.  I thought if I set:

<sld:VendorOption name="followLine">false</sld:VendorOption>

It would ignore this parameter.

How do I make my labels ignore followline?


texts: Difference Between Radiometric calibration and radiometric correction

texts:
I am new in remote sensing field.
I sometimes got confused when I read about image pre-processing.
Can someone please elaborate to me the difference and some example methods to do radiometric correction and radiometric calibration ?
Is converting DN values to Top Of Atmospheric (TOA) reflectance values from 2 sets of images (e.g Landsat 7) from different time counted as correcting or calibrating ?
Are DOS (Dark Object Subtraction) and Sun-angle correction parts of radiometric correction ?


texts: Layout expression for getting all features of layer within the current map view in QGIS 3?

texts:
In a map layout template, I want to insert an expression which gets all feature's names of a certain layer, which are contained in the current layout's view/extent. Is there a way to get an aggregate of all features of a layer which are contained in the current map view?

I've come so far, to check if a single feature of a layer is contained, but don't know how to go on from here:

contains(
         map_get(item_variables('Main Map'), 'map_extent'),
         nodes_to_points(geometry(get_feature_by_id('Corridors_86ef2ea9_d1a1_4d9f_9735_7b7b2fb54cb2', 2)))
)

I'm quite sure there is an expression combined with the aggregate function, which would enable me to concatenate all the names of the features in the current layout view and use this as a text item in my final layout.


texts:
I am using QGIS 2.14 and have problems with labeling a layer.

My layer is a MSSQL-Layer (from a GIS-Server) and the fieldtype I am using for labeling is a Sting (typename is nvarchar).If the string contains a forward-slash (/), then no label is displayed. Only objects without forward-slash in the labelfield are labeled.

I can use the same layer in QGIS 2.8 without problems on labeling (all object, including forwar-slash, are labeled).

For testing I created a local shp-layer with the fieldtype string (and typname string) and I have no problems on 2.14.

Have anybody an idea what´s going wrong on labeling my MSSQL-Layer.


texts: Manipulating Elevation Data from the USGS National Map

texts:
I am trying out @tylermorganwall Rayshader R library. Some instructions are given here by tylermw8 to get elevation data from the USGS National Map.

I am interested in obtaining data around the Golden Gate Bridge.
The data I download from from USGS National Map (n38w123) is, I believe, a one degree square tile which is 60 miles N/S and about 53 miles W/E. That's too much data for playing around with his rayshader library.

The bounding square of interest is:

(37.88, -122.56)   (37.88, -122.41)

(37.77, -122.56)   (37.77, -122.41)

How do I 'cut' this elevation data out of the large download from USGS? Is this where the library(sf) comes in handy?


texts: Checking if layer exists within particular dataframe using ArcPy?

texts:
I have a mxd with 3 dataframe: 'layers', 'geography' and 'book'.
'layers' dataframe has 2 layers: 'gridIndexFeatures' and 'zone'.
'geography' and 'book' datasets are empty.

I want to add the 'gridIndexFeatures' layer to 'book' dataframe by means of a  Python script, but I want to check if that layers was previously added so as to not add it again.

I wrote the following code:

def cargacapas(capa, alias, dataframe, mxd_path):
    arcpy.MakeFeatureLayer_management(capa, alias)
    DF = arcpy.mapping.ListDataFrames(mxd_path)[dataframe]
    layer = arcpy.mapping.Layer(alias)
    if not arcpy.Exists(layer):
        arcpy.mapping.AddLayer(DF, layer, "AUTO_ARRANGE")
    else:
        print "it exists so I don't add it again"
    mxd_path.save()
    del mxd_path

mxd = arcpy.mapping.MapDocument(r"U:\Proyectos\book1.mxd")
grid_path = r"U:\Proyectos\EDAR.gdb\gridIndexFeatures"
inFeatures = r"U:\Proyectos\EDAR.gdb\zone"
cargacapas(inFeatures, "Zona_de_ocupacion", 2, mxd)
cargacapas(grid_path, "Grid", 2, mxd)

My problem is that the EXISTS function detect 'gridIndexFeatures' exists in the 'layers' dataframe and return a True. But I need to detect if it exists in the 'book' dataframe.

Any help with this issue?

I have changed the code like this:

def cargacapas(capa, alias, dataframe, mxd_path):
    arcpy.MakeFeatureLayer_management(capa, alias)
    DF = arcpy.mapping.ListDataFrames(mxd_path)[dataframe]
    layer = arcpy.mapping.Layer(alias)
    if layer not in [lyr.name for lyr in arcpy.mapping.ListLayers(mxd_path,"", DF)]:
        arcpy.mapping.AddLayer(DF, layer, "AUTO_ARRANGE")
    else:
        print "it exists"
    mxd_path.save()
    del mxd_path

But it not seems to work propertly. Each time I run the code, it add the two layers, so it doesn't detect that those layers are already in the dataframe.


texts: KMZ file displays correctly in Google Earth but not inside Google Maps

texts:
I have a KMZ file that contains many placemarks inside.

When opened in Google earth the file displays the placemarks correctly. I can click on them and see the data for each placemark. See screen shot below which shows the placemarks appear on top:

I transferred the KMZ file to Google maps, and now I can't click on some of the placemarks. Only the top layer seems to work and I can see the placemarks behind everything but can't click them like I can in Google Earth.

I assume both programs function differently and that is why, but is there anything I can do to fix the KMZ file so it displays correctly inside Google Maps?


texts:
I'm working with a dataset with latitude, longitude and date-time. I have millions of points per day and I want to group them using a clustering algorithm.

I read this question Clustering geographic points but in my case I don't have an expected number of clusters. I explored HDBSCAN but it only accepts 2 features (lat lon in radians).

The idea is to get clusters of my datapoints regrouping nearest points in location and time together. So depending on the day, the output should change.

Is there a way to do a spatiotemporal clustering that includes the 3 features?


texts: How to use Select by Attributes query with semicolons?

texts:
I am trying to select multiple (50 or more) unique numeric records using ArcGIS 10's "Select by Attributes" using semicolons ';' after each record which it verifies as good to go.
Please see screen shot of where it verifies my expression "successful".

Yet when OK is hit it only selects the first record.

Do we really have to enter "xyz" = 1 OR "xyz" = 2 OR "xyz" = 3 etc. etc. to have multiple records selected?


texts: Snapshot live feed featurelayer by creating a copy of the data in ArcGIS

texts:
Let's say we have a feature layer on ArcGIS online that is constantly updated.  Examples would be the Recent Conditions in Air Quality (PM2.5), Satellite (VIIRS) Thermal Hotspots and Fire Activity, or in my particular case, the Recent Hurricanes, Cyclones and Typhoons layer.
These layers are constantly changing, but I want to create a snapshot of the data as is, so that I can have a record of say, Hurricane data from 2020, or Air Quality forecast from today.  I tried following the Copy and save layers tutorial for ArcGIS online, and I ended up with Hurricanes, Cyclones and Typhoons - Observed Position - 2020, but I quickly realized that this is just referencing the original data, it does not create a snapshot as a new data source.
I tried opening a projectin ArcGIS Pro and bringing in the layers, hoping that I could export the features to a new dataset and then publish that data as a new layer.  But when trying to go to export the data, the option is greyed out:

What can I do to export / copy this data and publish it as a new layer that does not update alongside its original source?
Surely there must be a way.  I have access to ArcGIS Online and Pro, and my eventual goal is to publish the layers on ArcGIS online and then consume them in an ArcGIS API for JavaScript app.


texts: How can I insert odd and even number in attribute table for selected feature in Arc Map?

texts: I work with addressing, in right site of street must be inserted even number, and in left site of street must be inserted odd number in attribute table. Is there any code to insert odd and even number spontaneously from field calculator?


texts:
I am using QGIS 2.12.3

I have two layers, one a complete grid of 1km squares, the other polygons representing urban land.

Whether I intersect clip or difference these two layers, some of the km squares are omitted from the new layer. It can't be a data issue as the missing squares are different depending on which operation I undertake, although they appear consistent for each operation.

Can anyone explain - I must admit I find the geoprocessing tools can often ignore a polygon on an layer with no obvious reason (yes I have tried validating the layers with no errors).

Note how the same polygon interacts properly with many other squares.

The polygon does fail validation but I have been around the periphery of failed squares and can see no faults.


texts:
I have some images at a local coordinate system with the following specifications:

TIFF header file: -gdalinfo im.tif

Driver: GTiff/GeoTIFF

Coordinate System is:

PROJCS["Eastern Ridge Project Grid GDA94",
GEOGCS["GDA94",
    DATUM["Geocentric_Datum_of_Australia_1994",
        SPHEROID["GRS 1980",6378137,298.2572221009113,
            AUTHORITY["EPSG","7019"]],
        AUTHORITY["EPSG","6283"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4283"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",120.95],
PARAMETER["scale_factor",0.9999],
PARAMETER["false_easting",446904.02],
PARAMETER["false_northing",2879827.84],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]

Origin = (326250.000000000000000,297750.000000000000000)

I want to be able to transform the images to UTM/EPSG:28350 using gdalwarp.
I use the following command for this purpose:

-gdalwarp -t_srs EPSG:28350 input.tif output.tif

and

-gdalwarp -s_srs '+proj=tmerc +lat_0=0 +lon_0=120.95 +k=0.9999 +x_0=4
46904.02 +y_0=2879827.84 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
-t_srs EPSG:28350 im.tif out.tif

ERROR 1: Translating source or target SRS failed:
'+proj=tmerc

When I open the image file in QGIS, it shows that projection is 28350 but coordinates aren't correct and show the previous system as if no transformation has been done.
I also put the information above in a text file and used gdalwarp but no luck.

I have two questions:

Are these parameters sufficient to define a projection? If yes, how can use it in gdalwarp or is there any other possible solution?


texts: I have added google map with PyQt5 and PyQGIS and I want to overlay points and images on map(specifically in python code) but I don't know how to do that. I am using PyQt5 and QGIS3.


texts:
I have file : test.kml and I try to open it with GDAL/OGR (1.6) librairie in Python (2.5).

But I don't succed can you help me ?

This is my code :

from osgeo import ogr

adresse = 'test.kml'
driver = ogr.GetDriverByName('KML')
datasource = driver.Open(adresse)
layer = datasource.GetLayer()
layerDefn = layer.GetLayerDefn()
featDefn = layer.GetLayerDefn()
feat = ogr.Feature(featDefn)
nbFeat = layer.GetFeatureCount()

print nbFeat

  ERROR Traceback (most recent call
  last):   File "C:\Documents and
  Settings\Guilhain\Mes documents\My
  Dropbox\dev\test.py", line 6, in

      layer = datasource.GetLayer() AttributeError: 'NoneType' object has
  no attribute 'GetLayer'


texts:
I am using GDAL 1.91 on a Windows 7 machine.  I have been mosaicing large numbers of rasters together and most of the time this process works exactly as expected.  However a very few of my mosaics look incomplete.  It is as though GDAL_merge simply ignored some of my input files (see image - notice the top lefthand corner):

To try and eliminate this problem I have:

checked that the missing rasters are valid and exist in my file system (correct on both counts - they display individually without error)
checked my input string and that definitely includes the missing rasters.
made all the input rasters the same type (e.g. Float32)
explicitly specified the no_data value to ignore in the inputs
EDIT: all the inputs are also the same resolution and running with verbose output gives no clues of any errors either

The mosaic process involves joining a maximum of 395 rasters each of which is 980KB and the resultant files should be about 390MB.  The process works very well most of the time, so this suggests that it is not a memory issue (I have 6GB of RAM and am using the 64 bit version of GDAL).

The error is not isolated to only one mosaic.  I have seen this elsewhere before and with GDAL 1.8 (which is why I upgraded to 1.91). In this particular case, the result using GDAL 1.8 was almost the exact opposite (i.e. the top left hand corner of this area did have data and the bottom left was empty - thus proving my data are valid but mystifying me further as to what is going wrong).

Changing the order of my input files does sometimes make a difference but usually does not entirely solve the issue.

Is this a bug (this is my gut feeling at the moment)?  Otherwise, does anybody have any suggestions for what might be going wrong and/or work arounds?

Many thanks!

UPDATE:
I have been studying my rasters in QGIS looking for differences and the ONLY difference I can see is that, although all the rasters were created in exactly the same way and are all the same resolution, same projection and all geotiffs, QGIS quotes the driver as VRT/Virtual Raster for the ones which gdal_merge is failing.  QGIS correctly quotes the driver as GeoTiff for the rest!?!  Gdalinfo on the other hand, quotes the driver as GTiff/GeoTiff for all the files.

So - there is SOMETHING different about the tiles but what and why is not obvious... more digging to follow...


texts:
I am trying to load a GeoJSON layer from my GeoServer WFS and have followed the many question responses but still cannot get it to work.

If i enter the url into a browser the GeoJSON appears with all the data, but the layer will not load in a leaflet map.

I also receive this error in the Google Chrome console:XMLHttpRequest cannot load http://localhost:8080/geoserver/IT_Demo/ows?service=WFS&version=1.0.0&reque…_Demo:states&outputFormat=application/json&format_options=callback:getJson. No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'file://' is therefore not allowed access.

Any ideas what may be wrong with my code?

<!DOCTYPE html>

  <title>A Leaflet map!</title>
  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css"/>


texts: Openlayers 5 scaleline has gives wrong results

texts:
I am using Openlayers 5 (modules from NPM) and I am having a problem with a scaleline in my map.

// from https://epsg.io/2264.js
proj4.defs('EPSG:2264', '+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.2192024384 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs');

register(proj4)

opts = {
  extent: extent,
  center: [
    (extent[0] + extent[2]) / 2,
    (extent[1] + extent[3]) / 2
  ],
  maxResolution: 1777.77777777778,
  minZoom: 2,
  maxZoom: 22,
  zoomFactor: 2,
  zoom: 2.5
}

// ---------------
// I get different results in scaleline depending
// upon which of these projection lines is used

opts.projection = 'EPSG:2264';
//opts.projection = new Projection({code: 'EPSG:2264'});

// -------------

let view = new View(opts);
let map = new Map({
  view: view,
  layers: []
});
map.addControl(new ScaleLine({
    units: 'us'
}));

I get different scaleline results depending upon which projection line above that I comment out.

When the map first appears, using the first projection line above gives a 2 mile scaleline. Using the second one gives a 5 mile scaleline. This is the same length scaleline and map extent but a different distance label on the scaleline.  The second seems to be the correct one.

Shouldn't these give the same results?

Am I doing something wrong or should I be filing an issue?


texts:
I'm trying to execute this (it's JavaScript but it doesn't matter):

queryTask.execute(query, function (results) {

and I get an error:

  Failed to execute query

These are the parameters:

var queryTask = new esri.tasks.QueryTask(path);
var query = new esri.tasks.Query();

To show how are defined what I get with them, console.log shows this:

console.log('query', query);

outFields: Array(1)
0: "*"
length: 1
__proto__: Array(0)
returnGeometry: false
spatialRelationship: "esriSpatialRelIntersects"
where: "codigo_pos = 8150215"

console.log('queryTask', queryTask);

gdbVersion: undefined
normalization: true
source: undefined
url: "https://arcgis.XXX.com/arcgis/rest/services/XXX/XXX/MapServer/4"
_countHandler: ƒ ()
_errorHandler: ƒ ()
_executeForIdsHandler: ƒ ()
_extentHandler: ƒ ()
_handler: ƒ ()
_relationshipQueryHandler: ƒ ()
_url: {path: "https://arcgis.XXX.com/arcgis/rest/serv…XXX/XXX/MapServer/4", query: null}
__proto__: Object

Any idea?


texts: TimeSeries and Export to CSV Google Earth Explorer

texts:
I am fairly new to coding.
I have the below code to generate monthly NDVI data for each of the 15 regions.  I can do them all individually, but it does take some time.
Are there any method of shortening l the process?
Also then I open each chart and download in CSV individually.
Also any code that can help to shorten this process?
IE I do the print chart 15 times.
var roi = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(18.445709161313030,-33.301100219893542), {name: '00'}),
  ee.Feature(ee.Geometry.Point(18.3283136995220101,-33.496446074764073), {name: '01'}),
  ee.Feature(ee.Geometry.Point(18.3283136995220101,-33.363695665206592), {name: '02'}),
  ee.Feature(ee.Geometry.Point(18.358508029205840,-33.477911798216702), {name: '03'}),
  ee.Feature(ee.Geometry.Point(18.555547014896824,-32.963288018732484), {name: '04'}),
  ee.Feature(ee.Geometry.Point(18.784147795713817,-33.421796302215846), {name: '05'}),
  ee.Feature(ee.Geometry.Point(18.797906906920584,-33.084616837828598), {name: '06'}),
  ee.Feature(ee.Geometry.Point(18.736492301824644,-33.379593144167053), {name: '07'}),
  ee.Feature(ee.Geometry.Point(18.608583989022776,-33.304010136472058), {name: '08'}),
  ee.Feature(ee.Geometry.Point(18.606336164825368,-33.520022728206399), {name: '09'}),
  ee.Feature(ee.Geometry.Point(18.919324659544490,-33.313052298558752), {name: '10'}),
  ee.Feature(ee.Geometry.Point(18.239021712114891,-33.161666295743785), {name: '11'}),
  ee.Feature(ee.Geometry.Point(18.298810314491927,-33.220923915024990), {name: '12'}),
  ee.Feature(ee.Geometry.Point(18.519748581700082,-33.348942987836971), {name: '13'}),
  ee.Feature(ee.Geometry.Point(18.784555700711259,-33.282413353935581), {name: '14'}),
]);
Map.addLayer(roi);

var startDate = ee.Date('1982-01-01'); // set start time for analysis
var endDate = ee.Date('2020-12-30'); // set end time for analysis

var nMonths = ee.Number(endDate.difference(startDate,'month')).round();

var ndvi = ee.ImageCollection('NOAA/CDR/AVHRR/NDVI/V5').select('NDVI')
                  .filterDate(startDate, endDate);

var byMonth = ee.ImageCollection(
    ee.List.sequence(0,nMonths).map(function (n) {
      var ini = startDate.advance(n,'month');
    var end = ini.advance(1,'month');
    return ndvi.filterDate(ini,end)
                .select(0).mean()
                .set('system:time_start', ini);
}));

// plot full time series
print(
  ui.Chart.image.series({
    imageCollection: byMonth,
    region: table0,
    reducer: ee.Reducer.mean(),
    scale: 1000
  }).setOptions({title: 'NDVI over time 0'})
);


texts: Calculating a new Latitude & Longitude based on percentage

texts:
I'm trying to work out the best way to do the following:

  Given coordinates in decimal degress (e.g. [0.0, 0.0]), calculate a new position to "pan" to, based on a percentage.

Let's say my chosen percentage to pan by is 1, how should I cater for the case when the starting point is [0, 0]? Performing simple arithmetic (i.e. multiplication) would not give a correct result (because 0 * x = 0).

For example, if I want to move from [0, 0] decimal degrees, 1% north, I would increase X by 90 * 0.01: i.e. latitude 0.9N.

Or to put it another way, what's the best way for me to calculate a new latitude and longitude by growing/shrinking another by a percentage, assuming within correct boundaries [-90,90] and [-180, 180]?


texts: QGIS join attributes by location python error

texts:
Relatively new to QGIS. Trying to add an attribute from a polygon to all the points in a shapefile (about 800 points) that are within that polygon. Used "join attributes by location", point layer as the target layer, polygon layer as the join layer. I ran it and after about 20% was processed, it stopped and I received the following error:

TypeError: QgsFeature.setGeometry(QgsGeometry): argument 1 has unexpected type 'NoneType'

Looking for some troubleshooting tips on how to identify what is tripping me up here.


texts: Converting L.GeoJSON Leaflet to different projection?

texts:
so I have some awesome basemap raster layers that use L.TileLayer to load, but the GeoJSON I have for my country polygons is all Web Mercator. My tile map is WGS84 or EPSG4326 (I think).

But I couldn't figure it out... I tried loading the "conversion methods" in Leaflet documentation, but couldn't figure out how to use it for pointToLayer() or coordsToLatLng() functions that are available as options on L.GeoJSON. I kept getting 'function not found' for static conversion methods. If they're static, they should be available but it means that I'm not doing this the proper way.

Using ES6 Leaflet:

let countryLayer = L.geoJSON(this.countriesGeoJson, {
          onEachFeature: function(feature, layer) {} });

Why is sample code so difficult to find for projecting to a different coordinate system? Obviously the "crs" option on Leaflet doesn't work because my GeoJSON is generated from the wrong map (was generated with Web Mercator projection) and there seems to be no easy way to convert.

I was trying to find a way without downloading a plugin like Proj4Leaflet ,but is that the only way??


texts:
I often get this error message when running pyqgis scripts on QGS project in EPSG:54003 or 54004 (I also noticed also some other SRSs in the past):

ERROR 6: EPSG PCS/GCS code 54003 not found in EPSG support files.  Is this a valid EPSG coordinate system?

Indeed, EPSG:54003 is implemented in QGIS (v2.18) by default and I can normally work in the QGIS desktop in this SRS.

Sometimes the script goes through, sometimes it stops.

What is the reason for the error message?
Where are the EPSG support files located?


texts: How to make a proper unification function

texts:
Well i make a map with OpenLayers 4 in draw a a brunch of polygons an try to union them if they intersect.

i made this code :

    function unification(){

 var parser = new jsts.io.OL3Parser();
 var w=0;
 var z=0;
 var compteur = Object.size(vector_layer.getSource().getFeatures());

for (w=0;w<compteur;w++){

    if (!vector_layer.getSource().getFeatures()[0]){
    }
    else{

  var featureA = vector_layer.getSource().getFeatures()[0];
  var jstsGeomA = parser.read(featureA.getGeometry());
  for(z=0;z<compteur;z++){

    if (!vector_layer.getSource().getFeatures()[z+1]){
    }
    else{

      var featureB = vector_layer.getSource().getFeatures()[z+1];
      var jstsGeomB = parser.read(featureB.getGeometry());

      if (jstsGeomA != jstsGeomB ){
      var intersect = jstsGeomA.intersects(jstsGeomB);
      if (intersect){

        var union = jstsGeomA.union(jstsGeomB);
        featureA.setGeometry(parser.write(union));
        featureB.setGeometry(parser.write(union));
        vector_layer.getSource().removeFeature(vector_layer.getSource().getFeatures()[z+1]);
        vector_layer.setVisible(false);
        vector_layer.setVisible(true);
        //compteur++;
      }
      }
    }
  }
}
}
}

but when i use my fonction on several polygons the 4th is delete simply without giving the union result in another.

any one got an idea?

ty

Finaly i founs an answere but it isn t optimize. If someone can help me it will be really hepfull.

function unification(){

  var parser = new jsts.io.OL3Parser();
  var a=0;
  var b=0;
  var compteur = Object.size(vector_layer.getSource().getFeatures());
  var featureA ;
  var jstsGeomA ;
  var featureB ;
  var jstsGeomB ;
  var intersection ;

  for (a=0;a<compteur;a++){
if(vector_layer.getSource().getFeatures()[a]){

  featureA = vector_layer.getSource().getFeatures()[a];
  jstsGeomA = parser.read(featureA.getGeometry());

  for(b=0;b<compteur;b++)
    if(vector_layer.getSource().getFeatures()[b] ){
      if(a!=b){
      featureB = vector_layer.getSource().getFeatures()[b];
      jstsGeomB = parser.read(featureB.getGeometry());
      intersection = jstsGeomA.intersects(jstsGeomB);

      if (intersection){

          var union = jstsGeomA.union(jstsGeomB);
          var unionFinale = union.holes= [];//WTF?WTF?WTF???
          featureA.setGeometry(parser.write(union));
          featureB.setGeometry(parser.write(union));
          vector_layer.getSource().removeFeature(featureB);
          a=0;
          compteur--;
          break;
      }
    }
  }
}
  }

}


texts:
I have a plugin that creates a new layer, which works fine.
When I edit this layer using the QGIS GUI, and done editing, I get the question of if I like to save the layer. If I hit 'yes' (and I want to save the changes), QGIS crashes (QGIS died on signal 11ptrace: Operation not permitted.)

I suspect this has someting to do when I create the layer:

layer = QgsVectorLayer("point", "test_points", "memory")
provider = layer.dataProvider()
provider.addAttributes([QgsField("ID", QVariant.Int),
        QgsField("TYPE", QVariant.String),
        QgsField("lat",  QVariant.Double),
        QgsField("lon", QVariant.Double)])

layer.updateFields()

QgsProject.instance().addMapLayer(layer)
layer.editingStarted.connect(self.ayer_edit)
layer.editingStopped.connect(self.layer_edit_end)

What do I miss?

EDIT:
I suspect this behavior has to do with my creation of layers. I tried to open and create a layer from a ESRI shape. The layers displays correctly, but when editing and saving the changes, QGis craches.


texts: Disable Token Authentication for REST Services in ArcGIS Server

texts:
For test purposes I have activated Token Authentication through my ArcGIS Server.

Now I want to get back to Windows Authentication and I just simply can't find anywhere clear instructions to disable Token Authentication from the REST Admin side of it? Nor even with ArcGIS community/blogs.


texts: Create 1km Grid in QGIS from Canada Topo Data

texts:
I've been struggling with this all day and can't seem to get anywhere, even with the help of other forum posts.  I'm pulling vector data from the canadian gis system, and loading it into QGIS.  That works all fine and dandy and I'm darn near enjoying playing around with the maps.  The problem I have is that I cannot seem to create an accurate 1km square grid.

Generally it seems that a interval of around 0.1 or less is the only way to even have gridlines show up on the screen.  I'm doing this both through the print composer and using the research tools -> vector grid methods.  Both work exactly the same.  From what I understand, living in Saskatchewan I need the UTM 13N projection, so that's what I've been using.

Anyone have any ideas what I might be doing wrong?


texts: Select by location tool Vs layer.getFeatures(QgsRectangle(xmin,ymin...)))

texts:
I can't make sense of the difference between:

QGIS Select by location tool (Vector->Research tools->Select by Location") which gives the correct result :7 records match the extent:

Creating a custmised pyqgis tool that uses the extent of the canvas to create a QgsRectangle and uses the layer it gives 66 matching records...why? United Arab Emirates is nowhere close!

Grab the extent from the input parameter (canvas extent) and force them to float, then:

vlayer.selectByRect(QgsRectangle(xmin,ymin,xmax,ymax))
features = [f for f in vlayer.getSelectedFeatures()]

model_feedback.pushInfo("\nMatched selected features from extent intersection: {}\n".format(len(features)))

The countries layer is in WGS84, the project is in WGS84, the returned extent minx,miny,maxx,maxy is correct for the area I have around Denmark

Any ideas what's wrong or is there a bug at play?


texts: Time series linear interpolation method for cloud not processing properly in my MODIS dataset

texts:
I am trying to perform an interpolation to my 3D array (sub-block of MODIS time series dataset with a shape of (46,960,960)) using numpy linear interpolation method to replace cloud pixels with interpolated values. The interpolation function is designed not to create a new array but simply replace the nan values with interpolated values. The function is shown below:

def interpolate(array):
    arrN=np.array(array, copy=False)
    ix=-np.isnan(arrN)
    idxTrue=ix.nonzero()[0]
    valueTrue=arrN[-np.isnan(arrN)]
    countFalse=np.isnan(arrN).ravel().nonzero()[0]
    arrN[np.isnan(arrN)]=np.interp(countFalse,idxTrue,valueTrue)
    # Got this algorithm from a discussion forum

arrVeg holds my 3D array as described above. The next process will simply loop to each index in my 3D array to access the time series values of each index. Then it will perform interpolation to replace nan with an interpolated value.

dim=arrVeg.shape
for r in range(0,dim[0]):
    for c in range(0,dim[1]):
        interpolate(arrVeg[:,r,c])

I checked it with a sample dataset (1D array) and the interpolation seems to be working fine (see below). It also works with a small sample of (3,3,3) matrix. However, if I ran it with my real dataset I noticed that all the values in my array that are assigned to np.nan were just converted to zero.

>>> File=pd.read_csv("J:\\MODIS\\Test\\Test_Alg2.csv")

>>> pixr=np.array(File.PIXR)
>>> array([0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 3, 3, 1, 3, 1, 3], dtype=int64)

>>> ndvi=np.array(File.NDVI)
>>> array([6748, 6920, 6698, 6631, 6466, 4810, 4161, 3508, 3113, 2837, 3607, 3519, 3519, 2007, 4125, 4933, 6832, 5236, 2998], dtype=int64)

>>> filter=np.where(pixr==3,np.nan,ndvi)
>>> array([ 6748.,  6920.,  6698.,  6631.,  6466.,  4810.,  4161.,  3508., 3113.,  2837.,  3607.,  3519.,  3519.,    nan,    nan,  4933., nan,  5236.,    nan])

>>> interpolate(filter)
array([ 6748.,6920.,  6698.,  6631., 6466., 4810., 4161., 3508., 3113., 2837., 3607., 3519., 3519., 3990.33333333, 4461.66666667, 4933. , 5084.5, 5236., 5236.])

I checked the output of my processed data and the graph look like the image below. The blue line supposed to be the interpolated line, and the red line is the original NDVI value. Most of the pixels that are flag as cloud and was assigned by np.nan were just converted to zero. I checked my script for two days now and I haven't figured out yet what is causing this unusual behavior.

This is the function I am using to create my 3D array.

def array_stack(ndvi,qflag,xrows1,xrows2,yrows1,yrows2): #this should be list
    index=np.arange(0,len(ndvi))
    arrstack=np.empty(shape=(46,abs(xrows2-xrows1),abs(yrows2-yrows1))) #dtype=float
    for nv,q,x in itertools.izip(ndvi,qflag,index):
        arrndvi=gdalnumeric.LoadFile(nv)[xrows1:xrows2,yrows1:yrows2]
        #print "%s, x1=%s, x2=%s, y1=%s, y2=%s" %(nv,xrows1,xrows2,yrows1,yrows2)
        arrqflag=gdalnumeric.LoadFile(q)[xrows1:xrows2,yrows1:yrows2]
        mask_clouds=np.where(arrqflag==3,np.nan,arrndvi)
        arrstack[x,:,:]=mask_clouds
    return arrstack

Snapshot of my processed image. If you look at the table on the right side the value zero are the cloud pixels. If you look on the image, the grayish color are the cloud pixels also.


texts:
What is the best way to cite QGIS in a formal publication?  I mean, what text should be included in the References / Literature Cited section of a paper or report.

I'm also specifically referring to the software itself, rather than the user's manual.  I've looked in the user's manual and in Help->About, but I can't find an answer.


texts: R: plotting 2 polygon shapefiles, problems with CRS & extent

texts:
I am trying to create a map using 2 shapefiles, both of which have been subset. I am plotting municipalities of a census district in Canada.

For the census district, I am using the [census division boundary file from StatCan] (http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2016-eng.cfm). For the municipalities, I am using [CanVec files from Natural Resources Canada] (http://geogratis.gc.ca/api/en/nrcan-rncan/ess-sst/23387971-b6d3-4ded-a40b-c8e832b4ea08.html#distribution).

I have subset both files to cover the census district of interest and want to plot the CanVec administrative features (towns) on top of the district.  Individually, both plot okay, but

plot(towns,
     add = TRUE)

does not add the towns.  Nothing happens, no error.

I checked the CRS and found that my 'cd' (Census Division) uses NAD83. Documentation for the CanVec file, subset to 'towns', says it also uses NAD83, but my output looks like this:

proj4string(towns)
"+proj=longlat +ellps=GRS80 +no_defs"

So, I ran:

proj4string(towns) <- CRS("+proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675
           +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000
           +datum=NAD83 +units=m +no_defs +ellps=GRS80
           +towgs84=0,0,0")

Now, CRS looks normal, but they still did not both plot, so I checked extent and found they don't overlap at all.

extent(cd)
extent(towns)

> extent(cd)
class       : Extent
xmin        : 5986825
xmax        : 6267206
ymin        : 1322528
ymax        : 1474248
> extent(towns)
class       : Extent
xmin        : -94.81746
xmax        : -91.4896
ymin        : 48.51538
ymax        : 49.36678

My search of SE got me this far, but I don't know where to go from here.


texts: How do I export labels as editable text objects to be used in Illustrator? QGIS 1.8


texts: Is there a way to get more than 5 digits of precision when using a double in a QGIS memory layer?

texts:
addressLayer = QgsVectorLayer("Point?crs=epsg:4326&field=address:string(50)&field=city:string(50)&field=state:string(50)&field=zipcode:string(50)&field=elevation:double(20,8)&field=latitude:double(20,8)&field=longitude:double(20,8)&index=yes", "Search Results", "memory")

This is the line I used to make the layer. I have tried other combinations such as double(15, 8) etc. but the layer still defaults to length of 20 and precision of 5 for the double data type. I can't seem to find any documentation on this behavior.

Any ideas?


texts: Dividing two or more rows in single attribute table?

texts:
I want to divide rows in my attribute table. I need to divide the lengths of rows having specific IDs. For instance the length of ID 21 divided by the length of ID 11. An this calculation execute for all rows in the table This can be easily solve with the field calculator. But what if I want to keep dividing all the rows (e.g. id 22/12, 23/13, 24/14, etc...)

Is there a way to automatically do this calculation for all rows?


texts: How to change page size of pagelayout for printing custom size?

texts:
I have a A3 printer (size: 297mm x 420mm). But I am wanting to print my pagelayout with paper size is: 265mm x 380mm

I tried to change the page size of my pagelayout to 265mm x 380mm by using the statement:

pageLayout.Page.PutCustomSize(26.5, 38);
 pageLayout.Page.Orientation = 2;

And, I tried to setup printer properties to A3 size or custom size (265mm x 380mm) but it is not successfull (it shows incorrect position).

Do you know how to fix it?


texts:
Seems to be a recurring nightmare here, but it seems ESRI UK gave us a bum steer.  You do not need to install ArcObjects .NET SDK (and, of course, the obligatory VS IDE), but you do need to have an installation of one of the desktop products in order to be able to use the dlls.

Heres a link to the help page, but it appears you need to have ArcMap/ArcCatalog/etc installed in order to fully be able to deploy any application you've developed using ArcObjects?  Is this insanity on my part?

I would have thought the easiest way would be to bundle the dlls, but am told you can reference the GAC dlls (Not too sure I understand this so a steer would be appreciated).  However, ESRI docs seem to clearly state they cannot be deployed without an appropriate ArcGIS Desktop application installed.

Surely this is just for custom Desktop tools and not AGS tools?

Any ideas?


texts: How to do spatial join in autocad without Map3d?

texts:
My colleague tells me (vaguely) that spatial join is possible in standard versions of autocad using reference manager and xref.  Does anyone know how to do this?  Or if not, how to do a spatial join of data in autocad.

(Map3d versions excluded.)

Edit: Yes, I mean the same kind of spatial join as in arcmap.  I will probably have some extended entity data attached to a bunch of lines.  I am wondering how you would join this to data from either the nearest, or all overlapping, points or polygons.


texts: How to perform weighted overlay analysis or fuzzy c means clustering for creating management zone delineation into agriculture Field?

texts: I have problem with creating management zone delineation map with multiple input variables such as NDVI map, NDRE map, soil map, elevation, slope, and aspects. I inserted into all variables in Google Earth Engine with ready unsupervised classification. So, my question is that I want to combine them and create zone delineation map. I found that in ArcGIS Pro they have fuzzy overlay and weighted overlay spatial analyst tools. So, do they any similar approach in the GEE?


texts:
I have a raster layer that covers a large area (three cities). When I load the raster layer into a 3D visualizer, there are multiple areas with a drastic change in the slope rather than a gentle smooth change in elevation.

Is there any way to smooth out the differences in the elevation or edit it ArcGIS?


texts: How to make survey lines within polygon?

texts:
I am trying to figure out how to make survey lines within a polygon.

So first I make a polygon. Now I need to fill this polygon with lines, where I can choose the spacing between the lines, and the angle of the lines. All the lines have the same angle and spacing.

I am trying to do it in QGIS, but I also have ArcGIS for Desktop.

Does anybody have an idea of how to do this?


texts: Select Layer by Location - Arcpy - Why are my inputs being interpreted as invalid?

texts:
I am working on a tool in a python toolbox and am totally confused about what is going on and causing my select by location operation to fail (Error 000368).

I am trying to select the features in one layer that are completely within a specific polygon.
Here is the line where the error is occuring:

for (index,item) in enumerate(TempIDs):
    lyrBoundaries.definitionQuery = "[" + parameters[1].valueAsText + "]" + " = " + "'" + str(item) + "'"

    for (index2,item2) in enumerate(Spec_Layers):
        currentLayer = arcpy.mapping.Layer(str(item2))

        arcpy.SelectLayerByLocation_management(currentLayer,"WITHIN",lyrBoundaries)

        arcpy.CalculateField_management(currentLayer,parameters[3].valueAsText,hopper,"VB","#")

The error is occurring in the arcpy.SelectLayerByLocation_management() line.  The text of the error message I'm getting is: (edited to remove the actual text of the path)

    Executing: SelectLayerByLocation "[FULL PATH OF currentLayer]" COMPLETELY_WITHIN GPL0 "0 Unknown" NEW_SELECTION
Start Time: Fri Mar 01 10:29:29 2013
Failed to execute. Parameters are not valid.
ERROR 000368: Invalid input data.

My question is really: Why is lyrBoundaries showing up as GPL0 and not as the path of the layer?  Beyond that, any tips or hints about why the inputs are being interpreted as invalid?


texts: How to read a feature class in an Esri Personal GeoDatabase using R?

texts:
I have feature classes in a personal geodatabase on a Win 8.1 machine (r 3.4.2, rgdal 1.2-15) and, as mentioned here, there is a requirement for Access users to also interact with them.

I have adjusted the code sample found here (which works perfectly for a file geodatabase):

require(rgdal)

# The input personal geodatabase
pgdb <- "C:/path/to/your/personalgeodatabase.mdb"

# Confirm the PGeo driver is available
subset(ogrDrivers(), grepl("PGeo", name))

Which returns:

   name                 long_name write  copy isVector
49 PGeo ESRI Personal GeoDatabase FALSE FALSE     TRUE

Then, at this point, ogrListLayers(dsn) fails:

# List all feature classes in a personal geodatabase
fc_list <- ogrListLayers(pgdb)

  Error in ogrListLayers(pgdb) : Cannot open data source

I have tried adding path.expand:

fc_list <- ogrListLayers(path.expand(pgdb))

And I have also confirmed that for gdal (and so rgdal), that

  Personal Geodatabases are accessed by passing the file name of the .mdb file to be accessed as the data source name. On Windows, no ODBC DSN is required.

I successfully accessed the feature classes using R-ArcGIS but that requires an ArcGIS licence, and so does not really solve the problem.

Suggestions?


texts:
I'm trying unsuccessfully to perform updates on multiple linked tables using updatecursors.

The following code works as I would expect:

trackRows = arcpy.SearchCursor("tracks","Key_='b2879561-d751-e111-9e92-0050568e18a0'")
for trackRow in trackRows:
    arcpy.management.MakeFeatureLayer("tracks","currentTrack","Key_='" + str(trackRow.Key_) + "'")
    grXingRows = arcpy.SearchCursor(GradeCrossings,"LinkKey='" + str(trackRow.Key_) + "'");
    for grxrow in grXingRows:
        inspectionRows = arcpy.SearchCursor(Inspections,"LinkKey='" + str(grxrow.getValue("Key_")) + "'")
        for insprow in inspectionRows:
            grxDefRows = arcpy.UpdateCursor(GradeCrossingDefects,"LinkKey='" + str(insprow.getValue("Key")) + "'")
            for grxdefrow in grxDefRows:
                grxdefrow.LocationVal = 7
                grxDefRows.updateRow(grxdefrow)
            del(grxDefRows)
        del(insprow)
    del(grXingRows)
    arcpy.management.Delete("currentTrack")

However, when I change the 4th line to use arcpy.UpdateCursor, the script still runs, but the data rows that are updated in the inner loop don't get updated in the DB.

Any ideas?

EDIT: Having tested out the suggestions about deleting cursors and rows more thoroughly, I now have the following:

trackRows = arcpy.SearchCursor("tracks","Key_='b2879561-d751-e111-9e92-0050568e18a0'")
for trackRow in trackRows:
    arcpy.management.MakeFeatureLayer("tracks","currentTrack","Key_='" + str(trackRow.Key_) + "'")
    grXingRows = arcpy.UpdateCursor(GradeCrossings,"LinkKey='" + str(trackRow.Key_) + "'");
    for grxrow in grXingRows:
        inspectionRows = arcpy.SearchCursor(Inspections,"LinkKey='" + str(grxrow.getValue("Key_")) + "'")
        for insprow in inspectionRows:
            grxDefRows = arcpy.UpdateCursor(GradeCrossingDefects,"LinkKey='" + str(insprow.getValue("Key")) + "'")
            for grxdefrow in grxDefRows:
                grxdefrow.LocationVal = 7
                grxDefRows.updateRow(grxdefrow)
                del(grxdefrow)
            del(grxDefRows)
            del(insprow)
        del(inspectionRows)
        del(grxrow)
    del(grXingRows)
    arcpy.management.Delete("currentTrack")
    del(trackRow)
del(trackRows)

Still no luck!  I tried deleting the rows within the loop, and tried deleting the rows outside of the loop.  Neither solved my problem.

EDIT:  Now I feel a little silly.  I had a map service running which was maintaining a lock on BOTH of the feature classes I was trying to update.  I used arcpy.TestSchemaLock to test, and found that it returned false for both feature classes regardless of whether the first was using a searchcursor or an updatecursor.  Stopping the service solved my issue.  However there is still some funny behavior here that I'd like to understand:

According to the documentation, if the updatecursor can't get an exclusive lock, then it should error out...but my script ran just fine, no errors...why?
Even given that not erroring were a valid behavior, I still should not have been successful with the inner update EVER when the service was running, but it worked fine as long as it was the only updatecursor I was using.  Why?

EDIT: After sleeping on it and trying again, it seems the success from 2 days ago was a fluke.  Both yesterday and today I have been unable to get the script to update the database correctly.  The script is unchanged, and I have checked to ensure that there are no schemalocks on these feature classes (TestSchemaLock returns true).  I'm not sure what stars had aligned which are now misaligned, but I still have no solution to the problem.


texts: How do I query for all polygons within a distance from a point?

texts:
Call me stupid, but somehow I couldn't figure out how to do the following.

I've got a testing shapefile with 18 features, which I successfully imported into GeoDjango (1.3, running on Spatialite). As an practical example, I've taken a centroid of one of them (area.mpoly.centroid) and tried to find all areas within a 20 km distance (touching would be enough, needn't cover complete area).

Area.objects.filter(mpoly__dwithin=(centroid, D(km=20)))

However, this only gives me an empty list. Increasing the number to 200 km (in theory, this should cover all of my dataset) also yields nothing.

Am I doing something wrong?


texts: Polygons not showing in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.2. I have used the following procedure to creat polygon:

'Add XY data' to create points
'XY to line' to create line (but I could not see the lines)
'feature to polygon' to create polygons (but could not see the polygons)

Can anyone assist me what I should do to see the polygons?

Here are examples of my points: (25,12.5), (-17.540,-10), (-17.54,10), (42.5,-12.5), (51.403,-34.791), (51.403,-10).


texts: Google Earth Engine: convert image collection to image

texts:
I want make Supervised Classification for all images in imagecollection. I make one image with all bands, but i can't make classification, console give error:

  merged.sampleRegions is not a function

Maybe there's another way how to convert image collection to image with bands or what am I doing wrong ?

//there's my imagecollection with param
var imgcollection = ee.ImageCollection('COPERNICUS/S2')
  .filterBounds(roi)
  .filterDate('2017-01-05', '2018-01-05')
  .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT', 20));

//there's i make function which convert imagecollection to one image with bands
var mergeBands = function(image, previous) {
  return ee.Image(previous).addBands(image, ['B8', 'B4', 'B3','B2']);
};
//there i get one images with bands:8,4,3,2
var merged = imgcollection.iterate(mergeBands, ee.Image([]));

var newft = ee.FeatureCollection('ft:1k5OiHRqWMQ61dDqzVFVadbEB9T4DWp0_ppYPAl7X');

//there's i trying to classificate
var training = merged.sampleRegions({
  collection: newft,
  properties: ['Number'],
  scale: 30,
});


texts:
I'm creating a GUI in Python for a Google Earth Engine Application and this is the result so far.

I can't figure out how to insert the map generated by Google Earth Engine into the grey area in the picture, that is a canva into the root window.

Is there a method to do this (and to do this only when the button is pressed)?


texts: Non-HTML circle text labels on Leaflet map

texts:
I need to save image of Leaflet map with polygon, multi-polygons, circles and text below every circle. I have geojson file with circle and polygon coordinates (lat and lng). I managed to get image of Leaflet map with circles and polygons with Leaflet-image library (https://github.com/mapbox/leaflet-image), but it can't rasterize html and labels are made in html.

Is it possible to create non-html text labels below every circle which would be rasterized with Leaflet-image? Is there a way to show some text labels which can be rasterized with some other web mapping libraries, like Mapbox or OpenLayers?


texts: How to overlap two layers of points with different size

texts:
I am new in QGIS. Here is my problem :

I have two layers of points. In these two layers, my points have different degrees or size :

For the layer 1, the point "a" has a size of 3, and the point "b" has a size of 8. The color of these two points is red.
For the layer 2, the point "a" has a size of 5, and the point "b" has a size of 2. The color of these two points is green.

Problem : When I overlap these two layers, then if I place the layer 1 first, we won't be able to the see point "b" of the layer 2. Conversely,  if I place the layer 2 first, we won't be able to the see point "a" of the layer 1.

To fix the problem, I tried to go on :  "Propriété de la couche" -> "Symbologie" -> "Rendu de couche" -> "Entité" and I choose "Multiplier", to make points translucent, but the results are not satisfying me.

Is there a way to overlap these layers so that we place the points whose size is the smallest, above the points of larger size, despite their belonging to two different layers?


texts: QGIS Change file browser/chooser default location

texts: Is there a way to change QGIS file browser default location? Currently it points to bin folder of QGIS installation C:\Program Files (x86)\QGIS Brighton\bin.


texts: How to combine Lat Long in CSV file with zone name in shapefile without any common identifier

texts:
I have a CSV file that has Latitude, & Longitudes of various locations in NYC. I have related shapefiles that have the various zones names (like Battery Park, CentralPark etc)

Problem statement: I need to add a new column to the CSV file that has the zone name corresponding to each Lat Long pair.

Note:

The CSV file doesn't have any identifier common with shapefiles.
I am a newbie to R. Following is the code I have tried and after this I have no clue how to go about.

library(sqldf)
library(ggplot2)
library(rgdal)
library(data.table)
library(xtable)
library(tcltk)
library(rgeos)

points<-as.data.table(read.csv.sql(file,sql="select * from file where type=2 ",header=TRUE,sep=","))

zones<-readOGR("nyc_zones.shp", layer="nyc_zones")
proj4string(zone)
coordinates(points)<-~longitude+latitude
proj4string(points)
proj4string(points)<-CRS("+proj=longlat +datum=NAD83")
points<-spTransform(points,CRS(proj4string(zones)))
identical(proj4string(points),proj4string(zone))


texts: Setting a default value in a QGIS field stored in geopackage

texts:
I try to create a vector layer with some field values which should use a QgsDefaultValue.
It works fine if I use an in-memory layer (useGeopackage=False) but fails (no default value) if I use a geopackage layer.
If I manually set the default value of the geopackage layer using the GUI it work fine however.
A minimal example looks like:
useGeopackage = True
vl = QgsVectorLayer("Point", 'name', "memory")
if useGeopackage:
    # save layer to geopackage
    path = 'test.gpkg'
    options = QgsVectorFileWriter.SaveVectorOptions()
    options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteFile
    options.layerName = 'layerName'
    QgsVectorFileWriter.writeAsVectorFormat(vl, path, options)
    del vl
    # load the newly created layer
    vl = QgsVectorLayer(path+'|layername='+'layerName', 'name', 'ogr')

pr = vl.dataProvider()

# add field
f1 = QgsField("active", QVariant.Bool)
f1.setDefaultValueDefinition(QgsDefaultValue('true'))
pr.addAttributes([f1])
vl.updateFields()
QgsProject.instance().addMapLayer(vl)


texts: Raster diff: how to check if images have identical values?

texts:
Is there a means to check to see if any 2 given raster layers have identical content?

We have a problem on our corporate shared storage volume: it's now so big that it takes over 3 days to conduct a full back up. Preliminary investigation reveals one of the biggest space consuming culprits are on/off rasters that really should be stored as 1-bit layers with CCITT compression.

This sample image is currently 2bit (so 3 possible values) and saved as LZW compressed tiff, 11 MB in file system. After converting to 1bit (so 2 possible values) and applying CCITT Group 4 compression we get it down to 1.3 MB, almost a full order of magnitude of savings.

(This is actually a very well behaved citizen, there are others stored as 32bit float!)

This is fantastic news! However there are almost 7,000 images to apply this too. It would be straightforward to write a script to compress them:

for old_img in [list of images]:
    convert_to_1bit_and_compress(old_img)
    remove(old_img)
    replace_with_new(old_img, new_img)

...but it's missing a vital test: is the newly compressed version content-identical?

  if raster_diff(old_img, new_img) == "Identical":
      remove(old_img)
      rename(new_img, old_img)

Is there a tool or method which can automatically (dis)prove the content of Image-A is value-identical to the content of Image-B?

I have access to ArcGIS 10.2 and QGIS, but am also open to most anything else than can obviate the need to inspect all these images manually to ensure correctness before overwriting. It would be horrible to mistakenly convert and overwrite an image that really did have more than on/off values in it. Most cost thousands of dollars to gather and generate.

update: The biggest offenders are 32bit floats that range up to 100,000px to a side, so ~30GB uncompressed.


texts: Designing User Interface in QT designer and connecting with QGIS

texts:
I am trying to input a excel file into QGIS from the QT designer User Interface.

I want to select the attribute list from the UI and display the thematic maps from the selected attribute list to the UI.

As i am new to this field i don't know how to establish connection between QT designer and QGIS to display the maps on UI.


texts:
I have created 250 square meter grids using the Grid Maker tool.
I have given each grid a unique name.
I want to find and store bounding co-ordinates (Latitude and Longitude) of each polygon in the Grid table along with the unique grid name.
In short, I want to extract values that are displayed when you select "Get Info" upon Right-Click of any polygon.
Please check screenshots for reference.
I'm using MapInfo 15.


texts: How can I output an arbitrary number of file urls from a Synchronous Geoprocessing Service?

texts:
Scenario:

Custom Export Service. The service creates zip files with shps in them and return the url to the file. I understand if you intend to to run this in synchronous mode, you must have a tool parameter with output type as "File" or the synchronous service will immediately delete the output zip file.

When testing, I noticed if my output parameter was instead set to "String" and the output content was a JSON object, async worked fine, but synchronous would not persist the output file on disk.

Question:

From a synchronous gp service, How can I return an arbitrary number of urls to output files on disk?

I'm wondering if this is impossible with current configuration options since output "File" parameters must be defined at publish time and cannot be altered at runtime.  I could return a json object with all the paths, but only any paths defined in the tool will remain on disk after the service's response.


texts:
I am using QGis version 2.2.0 and I set the snapping option to the vertex, with a tolerance of 1.5 meter. When I draw a new polygon next to an old one, it seems like the vertices are automatically anchored to each other, but actually when I zoom in to a scale of, let's say, 50:1, I can clearly see that the two vertices have not been anchored. Instead of one single vertex, I have then two of them.

It's very annoying because this automatic error can create topological problems of overlapping and gaps among polygons. In fact I detect these problems when I use the topology checker.

Did you encounter the same issue? How can I avoid it?

Thanks


texts: Custom encoding in qgis-web-client wsgi search

texts:
Firstly I would like to note that my qgis-webclient and the wsgi searching are all working correctly as they should. After creating afew datatables for test purposes the url:
http://localhost:1234/wsgi/search.wsgi?query=&searchtables=&query=11&cb=stcCallback1002
returned the correct values. Now my problem is that the encoding used in the wsgi script is using an python file somedisk\somepath\OSGeo4W\apps\Python27\Lib\json\encoding.py. With this file the text returned is encoded in UTF-8.

Now my question is how can I change this encoding to something different, in my case WIN1250 (ISO/IEC 8859-2).


texts: Can QGIS be used for commercial purposes? In particular, the output of maps/pictures for commercial use in reports.


texts: Using a dictionary to find unique values from one field, populate second field with consecutive integers

texts:
I have an attribute table with three fields: Id, Name, and Id_2 and 18,701 records. I would like to populate Id_2 with consecutive integers for each unique value in the Name field. Prior to running the script, Id_2 is empty. Below is an example of what I would like the final attribute table to look like.

Below is the python script that I have developed. Note that I have built a dictionary using the Name field as the keys and the Id field as the values. I don't need to use the values, however, to accomplish this task, just the key. My code produces consecutive integers, however, the sequence doesn't reset with each unique value in the Name field. Rather, my code populates Id_2 with consecutive integers for the entire table, through record 18,701 (Not what I want!).

# import required python libraries or modules
import arcpy
import traceback

# settings
arcpy.env.workspace = r'D:/test'
arcpy.env.overwriteOutput = True

# variables
fc = r'D:/test/colors.shp'
fields = ['Id', 'Name', 'Id_2']
for f in fields:
  print f, f in [field.name for field in arcpy.ListFields(fc)]

# build dictionary
rnDict = {}
cursor = arcpy.da.SearchCursor(fc, fields)
try:
    for row in cursor:
        id = row[0]
        rn = row[1]
        rnDict.setdefault(rn, []).append(id)
    print 'dictionary successfully built'
except:
    print 'An error occurred - part 2'
    traceback.print_exc()
del cursor

# update Id_2 field with sequential integers
startNumber = 1
cursor = arcpy.da.UpdateCursor(fc, fields)  # type: object
try:
    for row in cursor2:
        name = row[1]
        if name in rnDict:
            row[2] = startNumber
            startNumber = startNumber + 1
            cursor2.updateRow(row)
    print 'completed consecutive numbering'
except:
    print 'An error occurred during calculating field'
    traceback.print_exc()
del cursor
# del rnDict


texts: Python for finding the GPS coordinates of a point x-distance between 2 other GPS coordinates

texts: I have 2 pairs of GPS coordinates A = (44.27364400,-121.17116400) and B = (44.27357900, -121.17006800) and a distance of 10m. I'm trying to find the GPS coordinates of the point that's 10m from A toward B. I know I can use haversine to find the distance between A and B coutesy of: https://stackoverflow.com/a/4913653/5369777 is there someway I can add or subtract the distance 10m to return the degrees for the new point?


texts: Selecting Features in Visible Viewshed Areas using QGIS?

texts:
In QGIS 3.2 I have a viewshed raster and a shapefile with point features

I'd like to select the points that intersect the visible areas of my viewshed.

Is there a way to extract a vector of my visible areas so that I can then Select by Location?


texts: Rails - inserting GeoJSONs into PostGIS-adapted database

texts:
I'm using Rails' PostGIS adapter and RGeo-GeoJSON gems to store/serve some simple geometries in the Postgres database, but I'm having a little trouble with my database seed file... I'd like to be able to use GeoJSONs to seed my database, but seem to only be able to use WKT format.

For example, I'm able to seed the database successfully using the following syntax:

Model.create([{geometry: 'LINESTRING(-72.81164550781249 43.627896481020855, -72.90228271484375 43.386951440524854, -72.9324951171875 43.204107493733176)', description: 'test'}]),

but I'd much rather be able to simply insert a single GeoJSON into the geometry column, rather than have to break up its attributes across different columns and reassemble it later.

In the past, I've done this with PHP like so:

ST_SetSRID(ST_GeomFromGeoJSON('$geoJson'), 4326)

which passed the entire GeoJSON into one column, but it seems to me that the PostGIS adapter can only handle lat/longs in WKT format.  Is that the case, or is there a way to insert a pure GeoJSON?


texts: import *.osm.bz2 to postgresql using OSMOSIS causes errors

texts:
I am trying to import osm.bz2 flie from http://download.geofabrik.de/.
I think I am almost there though somehow my CMD returns error messages

I downloaded osmosis-latest file from http://bretth.dev.openstreetmap.org/osmosis-build/osmosis-latest.zip

what I followed was http://learnosm.org/en/osm-data/osmosis/ and system environment and everything else were set.

Yet, even the tutorial command.. which is extracting school POI and save it to  school.osm is not working in my PC. I have bhutan-latest.osm.bz2 file and *.osm isn't pbf so I switch --rbf to --rx (read xml)

this is the entire code I used

  osmosis --read-xml bhutan-latest.osm --nkv KeyValueList="amenity.school" --wx schools.osm

and the error messages are..

org.openstreetmap.osmosis.core.OsmosisRuntimeException: Argument keyValueListFile for task 2-nkv does not exist.
and more..

so from this concept
I tried to import the osm file to my postgis and this is my code

  osmosis \ --rx file=bhutan-latest.osm \ --write-apidb host="localhost" database="worldmap"port="5433" user="postgres" password="postgres"

and it give me this error code

  /mnt/c/Program Files (x86)/osmosis/bin/osmosis: 94: exec: java: Permission denied

by the way, I used both windows command and Ubuntu in windows 10

the reason why I am trying to use osmosis instead of osm2pgsql is I only need building footprint from OSM data. I believe osm2pgsql does not have extract query? If osm2pgsql can extract building footprint only, then I do not mind using osm2pgsql. osm2pgsql works fine in my Windows Ubuntu.

it takes so much time though, using windows10 ubuntu would be a lot easier for a person like me who does not have much knowledge of computer science.

you have to tune default.style.


texts: DEM - calculate shortest/easiest path over mountain range

texts:
In QGIS, I have a mountain range with a valley on either side (DEM raster file). Is there software or a plugin that allows me to calculate the easiest pathway over the mountain?
The answers at Cost Distance Analyst In QGIS? are from 3 years ago and involves skills I've not yet acquired.
Has nothing been developed in 3 years for us novices? Does QGIS3 not yet offer anything?