texts:
I'm working with pycharm 3.7 on win 8.1 with anaconda.
Im trying to import qgis with no success.
after reading this link I succeeded to import qgis when runing python from terminal but not from IDE.  It's happen after reading python-qgis.bat from terminal and then execute python from the terminal as well.

phthon-qgis.bat:

@echo off
call "%~dp0\o4w_env.bat"
call qt5_env.bat
call py3_env.bat
@echo off
path %OSGEO4W_ROOT%\apps\qgis\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES
rem Set VSI cache to be used as buffer, see #6448
set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%PYTHONPATH%
"%PYTHONHOME%\python" %*

I tried also to add this file path as an enviroment variable in pycharm but it change nothing.

how can I import qgis from my ide?


texts: Improving join between shape and view in ArcMap?

texts:
I have a problem with a simple join and I don't know what to do:

I have a shapefile with 400,000 parcels, in this file I have only one column with the parcel id with an index.
I have a Materialized View in Oracle with only 2,000 rows, this view has 4 columns, one of them is the parcel id also with an index.
In ArcCatalog I added a new database connection to get the Materialized View.
In ArcMap I added the shapefile and the materialized view with no problems. Displaying the data of each layer takes less than a second.
In ArcMap I create a join by parcel id keeping only matching records.

At this step, when I open the attribute table it takes up to a minute to get the data and drawing the layer also takes a long time. Because of this, I can't use this mxd project at all.

What am I doing wrong?


texts: ogrinfo: Unable to use sqlite dialect with -oo GEOM_POSSIBLE_NAMES open option

texts:
I am having problems using the SQLite SQL dialect to read a CSV file using ogrinfo, but only if I also specify a geometry column using the -oo GEOM_POSSIBLE_NAMES open option. Suppose file test.csv contains:

  id, geom
  1, 000000000140000000000000004010000000000000

Then this simple SQL works fine:

$ ogrinfo -GEOM=yes -sql 'select id from test'  test.csv
INFO: Open of `test.csv'
      using driver `CSV' successful.

Layer name: test
Geometry: None
Feature Count: 1
Layer SRS WKT:
(unknown)
id: String (0.0)
OGRFeature(test):1
  id (String) = 1

But now I wish to reference the geometry, so I use the -oo GEOM_POSSIBLE_NAMES open option:

$ ogrinfo -geom=YES  -sql 'select id from test' -oo GEOM_POSSIBLE_NAMES=geom test.csv
INFO: Open of `test.csv'
      using driver `CSV' successful.

Layer name: test
Geometry: Unknown (any)
Feature Count: 1
Extent: (2.000000, 4.000000) - (2.000000, 4.000000)
Layer SRS WKT:
(unknown)
Geometry Column = geom
id: String (0.0)
OGRFeature(test):1
  id (String) = 1
  POINT (2 4)

So far, so good.  Now suppose that I am doing something more complicated in my SQL (like CASTing or using a CASE), so I specify the SQLite dialect.  Now the wheels come flying off:

$ ogrinfo -geom=YES -dialect sqlite -sql 'select id from test' -oo GEOM_POSSIBLE_NAMES=geom test.csv
INFO: Open of `test.csv'
      using driver `CSV' successful.
ERROR 1: Cannot create virtual table for layer 'test' : CREATE VIRTUAL: invalid SQL statement : CREATE TABLE "test"("id" VARCHAR,"geom" VARCHAR,OGR_STYLE VARCHAR HIDDEN,"geom" BLOB_GEOMETRY_XY_0, OGR_NATIVE_DATA VARCHAR, OGR_NATIVE_MEDIA_TYPE VARCHAR)
ERROR 1: In ExecuteSQL(): sqlite3_prepare(select id from test):
  no such table: test

Notice that I can specify the sqlite dialect by itself just fine, so long as I leave off the "-oo GEOM_POSSIBLE_NAMES" open option:

/ # ogrinfo -GEOM=yes -dialect sqlite -sql 'select id from test'  test.csv
INFO: Open of `test.csv'
      using driver `CSV' successful.

Layer name: SELECT
Geometry: None
Feature Count: 1
Layer SRS WKT:
(unknown)
id: String (0.0)
OGRFeature(SELECT):0
  id (String) = 1

I noticed that when I try both the dialect and open options together, the error arises from GDAL's attempt to create and load a temporary SQLIte table with the contents of the CSV -- a table that is later accessed by the SQL.  The creation of this temporary table fails, possibly because the column name: "geom" is being specified twice.  That sort of makes sense as the open option causes GDAL to add the specify column as a geometry column, but it appears to have also been added as a regular column as well.  I've tried this on several different Linux flavors, and with varying GDAL versions.  Here's a specific crime scene:

/ # uname -a
Linux 1c20000629d8 4.9.0-6-amd64 #1 SMP Debian 4.9.88-1+deb9u1 (2018-05-07) x86_64 Linux
/ # gdal-config --dep-libs
-ljson-c -lsqlite3 -L/usr/lib -lpq -lz -lpthread -lm -lrt -ldl
/ # gdal-config --version
2.1.3
/ # ls -la /usr/lib/libsqlite3.so.0.8.6
-rwxr-xr-x    1 root     root        837024 May 20 21:06 /usr/lib/libsqlite3.so.0.8.6

I'm sure there other ways to manipulate CSV via GDAL with SQLite SQL, such as the OGR Virtual Format driver, but I'm really curious to learn what I'm doing wrong here.  Specifically, why is that specifying the SQLite dialect is incompatible with the use of -oo GEOM_POSSIBLE_NAMES?

Do I need to somehow exclude the geometry column as a "regular" column in the temporary table so that it is treated exclusively as a geometry column?


texts: Understanding DEM of difference results in ArcGIS Pro

texts:
I am looking at channel change in a river in the UK and am using three DTM rasters to observe change over time. One of the methods I am using is to simply minus one DTM from another.
I have got the results but I am unsure how to interpret them.
I believe the numbers are in mm.
EDIT:
So I have made sure that the DTMs are in the same co-ordinate system (British National Grid) and are alligned. Looking at the properties, it suggests that the scale is in metres? So what do the numbers represent e.g. scale +614 to +573168 for the 1998 DTM.
Even if the scale is incorrect and in fact is in millimetres, does this suggest areas that are eroded could have been so at 18,000mm (or 18 metres)? This seems too high as the river itself is only 0.5 metres deep.
(for the DTM of difference, the scale is between +18662 and -18492).
DTM 1998
DTM 2009
DTM 2009 - DTM 1998
Example of the difference


texts: Which tool/reference for address geocoding?

texts:
I have to geocode several corporate sites (a big excel file), close to Paris (France). I would like to geocode this file based on the addresses.

Could you tell me the best tool for this task (accuracy, performance, ease of use, etc.). On my station, I have  Qgis, postgis, mapinfo
Is it better to use a web tool or a desktop GIS?
I guess I must have a reference database for the streets network?  How do I choose it? And how do I get it?


texts: Drawing graphic layer on client side with ArcGIS API for JavaScript?

texts:
I'm using ArcGIS API for JavaScript 2.8.

I have some problems with drawing graphics point on client side. The problem is, when I draw points by first time by JS api everything is OK. But when points exists and I try to realy fast pan my map window from left to right (for example) my graphic layer can't finish drawing before one pan is ended. After few this kind of operations I can see, that my graphics points aren't in the geographic position as they should be when I started to pan. My question is how can I block again pan untill my graphic is drawn properly. I think it is the only way to prevent described situation.
I would like to mention, that I'm using jQuery for this application and I tried many many tricks like setting flag and setting enableMapNaviagtion() after document ready, I have tried to setTimeout, but for now this doesnt works.


texts: Error creating order with Python API: "\"products\" is not a mapping type

texts:
I am getting the following error when I try to place a Planet Labs order through the Python API as described at https://www.planet.com/docs/referencev0/orders/:

{"message": "Request validation failed", "errors": [{"field": "data", "message": "\"products\" is not a mapping type: Does not implement dict-like functionality."}]}

I am using the following Python code.

import requests
import os

SESSION = requests.Session()
SESSION.auth = (os.environ['PL_API_KEY'], '')

payload = {
  "type": "download",
  "data": {
    "products": [
      {
        "scene_name": "20150420_141621_081c",
        "product_type": "visual",
        "provider": "ortho"
      },
      {
        "scene_name": "20150420_141621_081c",
        "product_type": "analytic",
        "provider": "ortho"
      }
    ]
  }
}

r = SESSION.post("https://api.planet.com/v0/orders/", data=payload)
print(r.text)

What am I doing wrong?


texts:
I have a shapefile of lines and I would like to know an Arcmap tool to create a new shapefile  with orthogonal lines from the midpoint of that existing line shapes. I know that it is possible to do it manually, but I have so many lines so it would be an impossible way.

Does anybody know an automate way to do it?

Yes, perpendicular lines that cross the middlepoint of the existing lines. An example would be like this:


texts:
I have licensed ArcStorm from Esri but have no idea what the product is.

I have found no documentation about it from Esri.


texts: Changing style of points while moving them in QGIS

texts:
I have 1000 points in QGIS.
I want to move them one by one and I want to know which of those I have already moved.
I want to change style of point that I moved before.
How can I do that?


texts: arcpy.SelectLayerByAttribute_management fails on Feature Class

texts:
I'm trying to streamline contour production with a script tool. Typically we make the contours to bring into AutoCAD Map where we need the "Elevation" and "Layer" fields populated so they import correctly.

Basic steps in the tool:

-create a focal statistics raster

-contour the Focal Stats raster

-add CAD fields

-calculate the contour values to the "Elevation" field that Auto CAD recognized

-Calculate destination CAD Layers (intermediate or index)

I'm throwing an error when I try to do arcpy.SelectLayerByAttribute_management at the end of the script below.
the error states that "Parameters are not valid.The value cannot be a feature class".

OK. so if I can't use a Feature Class, what function should I use instead that can accept a FC input?

#Import stuff
import os
from arcpy import env
from arcpy.sa import *

#where am I saving?
path = (str(arcpy.GetParameterAsText(0)))

#What layer are we using ?
mxd = arcpy.mapping.MapDocument("CURRENT")
ras1= arcpy.GetParameterAsText(1)

#name of contour FC or shp
fc_out= arcpy.GetParameterAsText(2)

#zfactor
z=arcpy.GetParameterAsText(3)

#contour interval in units of z factor
interval=arcpy.GetParameterAsText(4)

#Focal Statistics Raster (hopefully this is just memory and goes away....confirm this later)
focalras=FocalStatistics(ras1, NbrCircle(3, "CELL"),"MEAN", "NODATA")

#do teh contour
Contour(focalras,str(path)+"\\"+str(fc_out),interval,0,z)

#Add CAD Fields

if ".gdb" in path:

    try:
        #Add CAD Fields
        arcpy.AddCADFields_conversion(str(path)+"\\"+str(fc_out),"ADD_ENTITY_PROPERTIES","ADD_LAYER_PROPERTIES","NO_TEXT_PROPERTIES","NO_DOCUMENT_PROPERTIES","NO_XDATA_PROPERTIES")

    except:
        # If an error occurs while running a tool print the message
        print arcpy.GetMessages()

#Calculate attributes in new contour FC - still working on this.....

    contours = str(path)+"\\"+str(fc_out)

    #copy "Contour" field to newly added "Elevation" CAD field
    arcpy.CalculateField_management(contours, "Elevation", '!Contour!', "PYTHON_9.3")

#--------------WORKS FINE UP TO HERE------------

    #Populate Indedx Contour layer name
    arcpy.SelectLayerByAttribute_management (contours, "NEW_SELECTION", "mod(Contour,10)=0")
    arcpy.CalculateField_management ( contours, "Layer", "\"ETOPO-IND\"", "PYTHON_9.3")

    #Populate Intermediate Contour layer name
    arcpy.SelectLayerByAttribute_management (contours, "NEW_SELECTION", "mod(Contour,10)>0")
    arcpy.CalculateField_management( contours , "Layer", "\"ETOPO-INT\"", "PYTHON_9.3")

else:
    print "SHP output - no CAD fields added"


texts: Change variables in QGIS 2.12 using python

texts: Now it's possible in QGIS 2.12 to add your own variables, I wonder how I can add new variables or change variables using the Python console.


texts: ListFeatureClasses overriding ListRasters in ArcMap script

texts:
What is preventing the script from identifying the spatial reference for the .tif files?  It appears that the script below lists the fc first, then gives the false spatial reference to the subsequent .tif files in the list.  Similarly, if there's a fc with an undefined spatial reference, the .tif files show "unknown".

arcpy.env.workspace = in_workspace
feature_classes = arcpy.ListFeatureClasses()
for fc in feature_classes:
    spatial_ref = arcpy.Describe(fc).spatialReference
    if spatial_ref.name == "Unknown":
        arcpy.AddError("{0} has an unknown spatial reference".format(fc))
    else:
        arcpy.AddMessage("{0} : {1}".format(fc, spatial_ref.name))
rasters = arcpy.ListRasters("*", "TIF")
for raster in rasters:
                arcpy.AddMessage("{0} : {1}".format(raster, spatial_ref.name))


texts: Conditionally combine two raster with raster calculator

texts:
I can’t figure out how to do this in QGIS. I tried gdal and SAGA raster calculator in QGIS. Here is the situation: I have two raster for the same area in the same projection with the same extent and the same resolution.

Raster A (has integer values from 1 to 5)

Raster B (has values from 1 to 7)

The calculation should do the following: if a cell in raster B is >= 3, make the cell for raster A + 1. So the result should be the whole raster A, but cell values in raster A where raster B is bigger than 3 should be +1.
My formula so far (for GDAL):

((A*(A>=3))+(B>=3))+(A*(A<3))

But the result is missing all values lower then 3. Anyone any suggestions?


texts: MassGIS Data Layer & PostGIS lat/lon not finding correct gemos

texts:
I am using the MassGIS House District data layer with srid 26986. I then query for the geom containing the lat/lon 42.21, -71.5 which should give me Boston but I got nothing. So I used OpenJump to visualize this and the point is rendering a few hundred miles south of the House District geoms:

I am guessing there is something wrong with the projection but according to multiple sources this is the correct srid: http://www.mass.gov/anf/research-and-tech/it-serv-and-support/application-serv/office-of-geographic-information-massgis/datalayers/overview.html.

Any pointers?


texts: Implementing RESTful Architecture on ZOO Project?

texts:
I read in the WPS 2.0 paper on Open GIS webpage that WPS 2.0 can be implemented with a RESTful architecture. I am working on a project using ZOO WPS.

Is there some way to incorporate the RESTful architecture in this?

According to my understanding so far, it seems that the framework is suitable only for SOAP architecture as all message formats are in xml.


texts: Is there an automatic way to connect attributes of a transformer to the attributes of a writer?

texts:
I am using FME 2015 in order to build a workflow as shown in the image:

As you can see in the figure I have an AttributeFilter which takes manually some value. I have to run this workflow multiple times for different values in the Attribute Filter. The issue I face is that once I modify the value (in AttributeFilter) then all the connections between the FeatureMerger and the writter (komvoi) are lost and I have to remake them by clicking..

My question is if there an automatic way to do this connection instead of manually.


texts: Clip to mask layer not adhering to expected width and height following CRM and resolution changes

texts:
I am attempting to take a simple opentopo heightmap and carry out the following in QGIS 3.4.1

- Reproject to UTM 31N
- Set cell size resolution to 5mtr
- Clip to a shapefile polygon feature exactly 40960 x 40960mtr

So why when using following command in Python console:

import os
os.system(r'''gdalwarp -t_srs EPSG:32631 -tr 5.0 5.0 -of GTiff -cutline E:\Users\rosso\Arma\QGIS\UTM31N\Test23-11-18\shape.shp -crop_to_cutline "E:/Users/rosso/Arma/Terrain Heightmaps/Heightmaps Opentopo/output_srtm.asc" E:/Users/rosso/Arma/QGIS/UTM31N/24-11-18/converted2.tif''')

Do I get width and height not matching exactly to my shapefile feature I'm clipping to - (40960 x 40960)  - which at 5mtr resolution size should be '8192' for both width and Height

If I break up the process into 3 stages:

Change CRM:

import os
os.system(r'''gdalwarp -t_srs EPSG:32631 -of GTiff "E:/Users/rosso/Arma/Terrain Heightmaps/Heightmaps Opentopo/output_srtm.asc" E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/crm.tif''')

Change Resolution:

import os
os.system(r'''gdalwarp -tr 5.0 5.0 -of GTiff "E:/Users/rosso/Arma/QGIS/UTM31N/25-11-18/crm.tif" E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/res.tif''')

Clip to Shapefile feature mask layer:

import os
os.system(r'''gdalwarp -of GTiff -cutline E:\Users\rosso\Arma\QGIS\UTM31N\Test23-11-18\shape.shp -crop_to_cutline "E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/res.tif" E:/Users/rosso/Arma/QGIS/UTM31N/26-11-18/cut.tif''')

The width and height adhere to what I expect:

I would prefer the process to be completed in one command, rather than splitting them out as I want to reduce the amount of processing on the original raster.


texts: OGR: Meaning of "2.5D extension as per 99-402"

texts:
I'm looking at the enum defs in ogr_core.h.  I see a bunch with comments that read:

  2.5D extension as per 99-402.

What in the world does that mean?  Presumably it's referencing a spec or RFC?

typedef enum
{
  ...
  wkbPoint25D = 0x80000001, /**< 2.5D extension as per 99-402 */
  wkbLineString25D = 0x80000002, /**< 2.5D extension as per 99-402 */
  wkbPolygon25D = 0x80000003, /**< 2.5D extension as per 99-402 */
  ...
}


texts: Is there an equivalent to ST_Rotate in spatialite?

texts:
I'd like to do some arbitrary rotation of spatial objects in spatialite (accessed through QGIS). In PostGIS I'd use something like:

...
SELECT geom AS ST_Rotate(geom, radians, originX, originY)

Is there an equivalent in spatialite?

Edit: after the answer and some more research, managed to find a way of mimicing ST_Rotate using ShiftCoords and RotateCoords, rotating by 45° around the centroid in the example below:

SELECT ShiftCoords(RotateCoords(ShiftCoords(GEOMETRY,
 -X(Centroid(GEOMETRY)), -Y(Centroid(GEOMETRY))), 45),
 X(Centroid(GEOMETRY)), Y(Centroid(GEOMETRY))) AS geom from TABLE


texts:
I have a problem here with installing ogr ..
I followed this very usefull tutorial for installing ogr on Windows systems: http://sandbox.idre.ucla.edu/sandbox/tutorials/installing-gdal-for-windows
It seemed to work, because I can call ogr from commandline afterwards.

C:\Users\root>gdalinfo --version
GDAL 2.1.0, released 2016/04/25

But when I try to use ogr2ogr it says that it can't find the drivers.
I'm on this for hours now but can't get it to work.
How do I register the drivers? The path is set correctly and I tried like thousands of different installations. All on a clean Virtual Machine which I reset after every failed installation ..

Please help me, I'm desperate :D

Edit:
This is the result I get when typing set g:

C:\Users\root>set G
GDAL_DATA=C:\Program Files (x86)\GDAL\gdal-data
GDAL_DRIVER_PATH=C:\Program Files (x86)\GDAL\gdalplugins

Which exactly matches my installation ..

Edit2:
When I'm trying to convert:

C:\Users\root>ogr2ogr -f "GeoJSON" "C:\Users\root\Desktop\Plugin\Shapefiles_Karte\Bundeslaender\Bundeslaender.geojson" "C:\Users\root\Desktop\Plugin\Shapefiles_Karte\Bundeslaender\Bundeslaender.shp" -overwrite
FAILURE:
Unable to open datasource `C:\Users\root\Desktop\Plugin\Shapefiles_Karte\Bundeslaender\Bundeslaender.shp' with the following drivers.
  -> netCDF
  -> PCIDSK
  -> JP2OpenJPEG
  -> DB2ODBC
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> OGR_SDTS
  -> S57
  -> DGN
  -> OGR_VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> NAS
  -> GML
  -> GPX
  -> LIBKML
  -> KML
  -> GeoJSON
  -> Interlis 1
  -> Interlis 2
  -> OGR_GMT
  -> GPKG
  -> SQLite
  -> ODBC
  -> WAsP
  -> PGeo
  -> MSSQLSpatial
  -> PostgreSQL
  -> MySQL
  -> OpenFileGDB
  -> XPlane
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDUMP
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> OGR_PDS
  -> WFS
  -> HTF
  -> AeronavFAA
  -> Geomedia
  -> EDIGEO
  -> GFT
  -> SVG
  -> CouchDB
  -> Cloudant
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> ODS
  -> XLSX
  -> ElasticSearch
  -> Walk
  -> CartoDB
  -> SXF
  -> Selafin
  -> JML
  -> PLSCENES
  -> CSW
  -> VDV
  -> TIGER
  -> AVCBin
  -> AVCE00
  -> HTTP


texts:
I am trying to make an application which will test if any of the segments (multi line layer acquired from GeoServer) intersects with polygon acquired using OpenLayers.Protocol.WFS.fromWMSLayer, and if it does add it to the map. After the polygon has been selected it's added to an empty vector layer which is then used as an input in WPS.Client.

The code looks like this:

int = client.getProcess('server', 'JTS:intersection');
   for (var i=0; i<v_layer.features.length+1; i++){
   int.execute({
       inputs: {
           a: selected.features[0].geometry,
           b: segment_layer.features[i].geometry
       },
       success: function(outputs) {
           map.layers[3].addFeatures(outputs.result);
       }
   });
}

This works, when I tested it on just a few segments, but since this layer contains about 200 of them it crashes the browser. I couldn't find any other WPS process that I could use other than JTS:intersection.

Any suggestions?

Thank you for your time.


texts: How to lock a shapefile in QGIS so only I can edit

texts: I am doing a project using shapefiles in QGIS. These will be sent to somebody who will view in QGIS and I want to put a lock on the data so that it can be viewed but not edited by somebody else without my consent. Can this be done


texts: Parametrized SQL Views with GeoServer and WFS

texts:
I ran into this problem that is (almost) nowhere on the internet and thought I'd share the solution. GeoServer supports parametrized SQL views for Web Feature Services (WFS) but the documentation leaves out a crucial fact (in fact its clear they cut and pasted the docs for the same page for WMS), since in a WFS all the requests are GET requests, they tell you to simply tack on to the query string a key value string of the form:

viewparams=keyOne:valueOne;keyTwo:valueTwo;

However in the GetFeature request, we have a POST request with an XML based payload and no obvious way to include the viewparams string. How does one do that?


texts: ArcObjects DeleteLayer leaving layer artifacts

texts:
I'm working on some code for a button in ArcGIS 9.2 using Visual Studio 2005. The process adds a point layer, performs a select by location, exports the attributes of the selected points, then removes the point layer. The code I'm using to remove the point layer at the end of the process is:

Dim pDL As IDataLayer2 = pPtLayer 'QI
pDL.Disconnect()
pMxDoc.FocusMap.DeleteLayer(pPtLayer)
pMxDoc.UpdateContents()
pMxDoc.ActiveView.Refresh()

But there are remnants of the point layer remaining after the process has completed. In the command window, a list of the layers indicates it is still being referenced. This is causing to button to perform incorrectly if it is used multiple times in a session. How can I remove the layer without leaving these 'artifact' references?

(Also posted at http://forums.arcgis.com/threads/54771-Data-References-Remain-After-Calling-DeleteLayer)


texts:
I need to write a code to read a polygon shapefile that contain three classes of land use (trees, urban, and crops) and I need to select only polygons of Trees.

Here is my code to apply Select by polygon that is not working:

class btn_read_LandUse(object):
    """Implementation for Addins_addin.button_LU (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        LU_layer = pythonaddins.OpenDialog("Please Select the land use file" ,False,r'','Select')
        desc = arcpy.Describe(LU_layer)
        layerName = desc.baseName
        arcpy.SelectLayerByAttribute_management(in_layer_or_view="LayerName", selection_type="NEW_SELECTION", where_clause=""""LU" = 'Crops'""")

The error:

Traceback (most recent call last):
  File "C:\Users\eng.adel\AppData\Local\ESRI\Desktop10.4\AssemblyCache\{FEC573D9-E048-4AEE-A331-8380627D272B}\Addins_addin.py", line 96, in onClick
    arcpy.SelectLayerByAttribute_management(in_layer_or_view="LayerName", selection_type="NEW_SELECTION", where_clause=""""LU" = 'Crops'""")
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 7182, in SelectLayerByAttribute
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Layer Name or Table View: Dataset LayerName does not exist or is not supported
Failed to execute (SelectLayerByAttribute).

Any help?


texts:
s-g filter had operated by GEE - Using Savitzky-Golay filter over an Image Collection, and I want calculate annual maximum value of image collection after s-g filter.
but there are not exist any bands in my collection of annual maximum value even through timestamp has added in image collection after s-g filter operation. I don't understand there's what problem in my code.
https://code.earthengine.google.com/d2e136d5ca0431d1fde2cd96794963f4
    // Run the SG solver over the series, and return the smoothed image version

var sg_series = runLength.map(function(i) {
  var ref = ee.Image(modis_res.get(ee.Number(i).add(half_window)))
  return getLocalFit(i).multiply(ref.select(indepSelectors)).reduce(ee.Reducer.sum()).copyProperties(ref)
})

    // This function adds a band representing the image timestamp.

var addTime = function(sg_series) {
  return sg_series.addBands(sg_series.metadata('system:time_start'))
                  .addBands(sg_series.metadata('system:time_end'));
};

var NDVI_sg = ee.ImageCollection.fromImages(sg_series)
print(NDVI_sg,"NDVI_sg")

var start_year = 2000;
var end_year = 2018;
var yearList = ee.List.sequence(start_year, end_year);
var yearImgList = yearList.map(function(year) {
  year = ee.Number(year);
  var ndvitempCol = NDVI_sg.filter(ee.Filter.calendarRange(year, year, "year"))
                  .filterBounds(roi);
  var NDVImax = ndvitempCol.max();
  NDVImax = NDVImax.set("year", year);
  NDVImax = NDVImax.set("system:index", ee.String(year.toInt()));
  return NDVImax;
});
var yearImgColNDVI = ee.ImageCollection.fromImages(yearImgList);
print("NDVIcollection", yearImgColNDVI);

var NDVI2000q = yearImgColNDVI.filter(ee.Filter.equals('system:index', '2000'));
var NDVI2000 = NDVI2000q.median()
print(NDVI2000,"NDVI2000")


texts:
I want to tile one large las tile into a few minor. But I just can't get it done.

My original las file have around 30 million points and I want to get 6 smaller las files with 5 million points in each.

What is that I don't know what exactly to put into "tile size" line.


texts: Do the elements in Print Composer have any identifiers in QGIS 2.18.14?

texts:
In QGIS 2:18:14, it is possible to refer to elements in the Print Composer?

How do I name and identify an element in another element?

E.g., on a label, call the contents of another label and use its content to do some calculations?

Let say: in a label, Label 1 and Label 2,  or Label 1 + Label 2, ...


texts: HTML5 local storage for tile maps and feature layers using ArcGIS Java Script API

texts:
I am developing a mobile arcgis application I want to store data locally so that when a user is offline he/she can use the  maps. I want to store tile layers,feature layers and a table data.

Any sample code appreciated

Platform : ArcGIS java script API and HTML5


texts: How to create a map with OSM and webatlas.de's WMS layer?

texts:
I spent hours trying to figure out the following problem, I hope someone can help me.

You can see here that this WMS has different Projections available:
http://www.webatlasde.de/arcgis/services/Maps4BW/MapServer/WMSServer?SERVICE=WMS&request=GetCapabilities

(here you can find even more information and an example map http://www.webatlasde.de/ArcGIS/rest/services/Maps4BW/MapServer )

How do I have to set up my map-object or layers so that the webatlas.de and OSM (or Google) are displayed in one Map with the correct projection?

Can't I use EPSG:3857 as an substitut for EPSG:900913?

In fact I can only get webaltas.de working without spatial distortion with the following configuration:

Proj4js.defs["EPSG:31467"] = "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";

var extent1 = new OpenLayers.Bounds(3364573, 5257369, 3635262, 5536211);

function Map_main(){

    map = new OpenLayers.Map("mapdiv",{

        projection: new OpenLayers.Projection('EPSG:31467'),
        displayProjection: new OpenLayers.Projection('EPSG:31467'),
        controls: [new OpenLayers.Control.Navigation(),
              new OpenLayers.Control.PanZoomBar(),
              new OpenLayers.Control.ScaleLine({ bottomOutUnits:'',bottomInUnits:'',maxWidth:150}),
              new OpenLayers.Control.MousePosition(),
              new OpenLayers.Control.LayerSwitcher()],
        units: "m",
        maxExtent: extent1
    });

    webatlasde = new OpenLayers.Layer.WMS("Webatlas", "http://www.webatlasde.de/arcgis/services/Maps4BW/MapServer/WMSServer", {
        layers: 'Halbtonvariante',
        format: 'image/png'
    });

    map.addLayer(webatlasde);

}

If I try this, the map is distorted

Proj4js.defs["EPSG:31467"] = "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs";

var leftBottom = new OpenLayers.LonLat(3364573, 5257369).transform(new OpenLayers.Projection("EPSG:31467"), new OpenLayers.Projection("EPSG:4326") );
var rightTop = new OpenLayers.LonLat(3635262, 5536211).transform(new OpenLayers.Projection("EPSG:31467"), new OpenLayers.Projection("EPSG:4326") );
var extent2 = new OpenLayers.Bounds(leftBottom.lon, leftBottom.lat, rightTop.lon, rightTop.lat);

function Map_main2(){

    map = new OpenLayers.Map("mapdiv",{

        projection: new OpenLayers.Projection('EPSG:4326'),
        displayProjection: new OpenLayers.Projection('EPSG:4326'),
        controls: [new OpenLayers.Control.Navigation(),
              new OpenLayers.Control.PanZoomBar(),
              new OpenLayers.Control.ScaleLine({ bottomOutUnits:'',bottomInUnits:'',maxWidth:150}),
              new OpenLayers.Control.MousePosition(),
              new OpenLayers.Control.LayerSwitcher()],
        units: "m",
        maxExtent: extent2
    });

    webatlasde = new OpenLayers.Layer.WMS("Webatlas", "http://www.webatlasde.de/arcgis/services/Maps4BW/MapServer/WMSServer", {
        layers: 'Halbtonvariante',
        projection: new OpenLayers.Projection('EPSG:31467'),
        displayProjection: new OpenLayers.Projection('EPSG:4326'),
        format: 'image/png'
    });

    map.addLayer(webatlasde);
}

Can someone tell me how to correctly change the projection for the webatlas.de WMS? Or do I even have to use its REST-Interface with OpenLayers (I tried but miserably failed)?

I hope I could make explain myself.

Have a nice day!
All the best
Gregor


texts: Cell Statistics standard deviation values incorrect

texts:
I'm having some issues with the Cell Statistics tool in a python script I am writing. I am using 30 years worth of global cumulative annual rainfall data as input and am getting strange values in some areas of my output raster.

In the image below the very bright areas have standard deviation values of around 2000. This is much higher than it should be. Besides these very bright areas the rest of the values in the raster look correct (tested by manually calculating standard deviation values). This is the standard deviation as calculated by my script.

Using the same set of input I calculated the standard deviation in ArcMap 10.2.2 using the cell statistics geoprocessing tool and got different results. The standard deviation seems to be correctly calculated across the entire extent of the raster.

I am providing statistics on the two raster in case it is helpful to anyone.

script raster:

min = 0
max = 6697.9404296875
mean = 75.83343676729601
std. dev. = 333.72000491046

arcmap raster:

min = 0
max = 1507.4833984375
mean = 35.804817248754
std. dev. = 84.35107381816201

There are no no data values under the issue areas that would be causing this. Why are only some of the values being calculated correctly when I process using arcpy? Any ideas about what is going on?

Edit: Added method from script where standard deviation is being calculated.

def calc_cv_raster(self, data_lst):
    print "\nCalculating coefficient of variation raster."
    std_dev = sa.CellStatistics(data_lst, 'STD', 'DATA')
    avg = sa.CellStatistics(data_lst, 'MEAN', 'DATA')
    cv = (std_dev / avg) * 100
    print "Calculated coefficient of variation raster."
    return cv


texts:
I'm trying to run a ogr2ogr command  with cmd. When I run the cmd the Osgeow command prompt just flashes and nothing happens. When I enter my ogr2ogr script directly to the Command Prompt everythinh works fine. My cmd script looks like this:

@ECHO OFF
CALL C:\OSGeo4W64\OSGeo4W.bat ogr2ogr -sql "SELECT * FROM 'D:\Testi2\own.SHP'.own WHERE Owner like 'Harr%'" D:\Testi2\own.shp D:\Testi2\All.SHP

What could be the problem?


texts: Error: Could not create and deploy an application with boundless sdk opengeo suite

texts:
I work on a web application with Boundless SDK opengeo 4.1,  I used the tutorial and followed the exchange and advice on this forum but I have some problems about adding plugins and deploying the apps. I need help if someone here can help me.

1-  I’ve already added all plugins I need for my project but I can’t add some plugins like gxp.Mesure, gxp.GoogleEarth 3D, gxp.FeatureGrid or Query Form which aren't enabled on the application. You can see below the application:

2-  I can’t make the package app for creating file.war to deploy with this application, I have a message error that you can see below:

3-  Is there a plugin that allows to export data like an excel format for the application with SDK clients?

Thanks !


texts:
I am using QGIS Version 3.16.4-Hannover and  Windows 10 but I also tested on Linux Mint 20.
I have a problem with several imported shapefiles from different sources. I do not know if I do something wrong or if it is simply not possible.
I wanted to add vector data from other sources to an existing QGIS-project that bases on a OpenStreetMap-XYZ-Layer (https://tile.openstreetmap.org/{z}/{x}/{y}.png). My project CRS is EPSG:3857 (WGS 84 / Pseudo-Mercator).
You can reproduce the problem by copying this sample of vector files from Euratlas: https://www.euratlas.net/tele/1700_vector_data_sample.zip
The file includes shape-files covering Switzerland/Savoy.

I added, for example, seas.shp (from the above mentioned ZIP) to my project. The file's CRS is strangely recognized as "ESRI:54004 - World_Mercator". After importing, the layer is misaligned, its offset is about 4 kilometres to north-east. See screenshot, blue layer.

I tried to correct this by manually switching the layers CRS to EPSG:3857. In result, the misalignment grew! Now, the layer's offset is about 20 kilometres to south-east. See screenshot, orange layer.

Here's the original seas.prj:
PROJCS["World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",0],UNIT["Meter",1]]

What am I doing wrong?
Is there any way to get the vector layer aligned to my project in EPSG:3857, thus to the OpenStreetMap-Layer?
Or do I misunderstand something important?
Is it simply not possible?
If so, could you explain why?


texts: ogr2ogr merge two GeoJSON to one GeoJSON

texts:
Is ogr2ogr capable of merging two GeoJSON files to a new one?
It is possible to merge shapefiles to one GeoJSON using the VRT driver. But when trying to do the same using GeoJSON files ogr2ogr fails.

VRT is:

<OGRVRTDataSource>
<OGRVRTUnionLayer name="union">
<OGRVRTLayer name="Bundeslaender">
<SrcDataSource>C:\Users\root\AppData\Local\OGR_Contextmenue\Bundeslaender.geojson</SrcDataSource>
</OGRVRTLayer>
<OGRVRTLayer name="standorte">
<SrcDataSource>C:\Users\root\AppData\Local\OGR_Contextmenue\standorte.geojson</SrcDataSource>
</OGRVRTLayer>
<GeometryType>wkbGeometryCollection</GeometryType>
</OGRVRTUnionLayer>
</OGRVRTDataSource>

And the ogr command:

ogr2ogr.exe -f "GeoJSON" "C:\Users\root\Desktop\Plugin_new\Shapefiles_Karte\Bundeslaender\Bundeslaender.geojson" "C:\Users\root\AppData\Local\OGR_Contextmenue\vrt.vrt" -overwrite

The error ogr2ogr gives me is:

ERROR 1: Failed to find layer 'Bundeslaender' on datasource 'C:\Users\root\AppData\Local\OGR_Contextmenue\Bundeslaender.geojson'.
ERROR 1: Failed to find layer 'standorte' on datasource 'C:\Users\root\AppData\Local\OGR_Contextmenue\standorte.geojson'.

To clear out: the files are existing at the path given in the vrt. I double checked it. Also the path where the merged GeoJSON file should be created exists and is accessible.
Doing the same with Shapefiles is no problem!
So what am I missing here?


texts: Changing fields name of an attribute table in a layout with PyQGIS 3

texts:
I have an attribute table in my layout, defined as follow :
layout = QgsPrintLayout(project)
layout.initializeDefaults()
QgsProject.instance().layoutManager().addLayout(layout)

pdf = QgsProject.instance().mapLayersByName('plan_de_ferme')[0]
pdf_table = QgsLayoutItemAttributeTable.create(layout)
pdf_table.setVectorLayer(pdf)
pdf_fields = ['champ','sup','supepnd']
pdf_table.setDisplayedFields(pdf_fields, True)

I would like to rename the fields name displayed in the layout without changing the fields name in the layer.
I tried :
pdf_headers = ['Champ','Superficie (ha)','Superficie épandable (ha)']
pdf_table.setHeading(pdf_headers)

Like it's the way to do with a QgsLayoutItemTextTable, but it's not working because "AttributeError: 'QgsLayoutItemAttributeTable' object has no attribute 'setHeading'".


texts:
When I see some MapServer I usually get deeper and checking how it was done. Recently I found on one server where is avalanche map strange tms url. I haven't seen it before but probably some of you guys will tell me what is this type. I know that 14 is zoom and rest is unknown for me. So if you know what is this tms url type, i will be thankful.

Before visiting this link you have to visit http://mapy.hiking.sk/?zoom=14&lat=6305837.94292&lon=2194989.75612&layers=00BFTTTTTTFFFFFTT to get cookies otherwise, this url below will not load.

https://mapy.hiking.sk/layers/laviny/14/000/009/088/000/010/772.png


texts: Splitting lines: Options for managing IDs and related records

texts:
I have roads in a feature class (Oracle 10.3 geodatabase).

Roads are split into separate lines where they intersect. This is an established business rule that cannot be changed.

Each road has related work order records/points (in a separate system).

Scenario:

A new road (#200) has been built that intersects an existing road (#100).

Road #100 needs to be split into two separate roads at the intersection of road #200 (using ArcMap).

Question:

What are the options for managing the the road IDs and the related work order records?

(I'm thinking in a generic/business-rule sense; I'm not looking for technical ArcMap instructions.)


texts: Why does On-the-fly reprojection not line up Google layers in UTM?

texts: My UTM shapefile lines up with Google layers if the project CRS is set to WGS84 mercator projection (including default 3857 and Google).  But if I set the project CRS to UTM Zone11, they are slightly shifted.  I've tried loading them and setting the project CRS in various orders - nothing worked.


texts:
I created two presets on qgis and they work very well when switching from one to another, yet on qfield and when I choose a map theme some layers are not displayed especially those that have a ruled style.

I'm blocked at this step.

How do I fix this?


texts: Is QGIS aspect analysis description mistaken?

texts:
I've been checking the Aspect analysis tool and several sources indicate that aspect is measured from 0 to 360, starting North and going counter-clockwise. Nevertheless, I've been testing it on several regions, some of which I know quite well and the measurement is actually being done clockwise.

Am I right? Are the descriptions and tutorials that state the direction as being counterclockwise mistaken?


texts: Insert new features into Oracle DB using PyQGIS

texts:
I need to insert new objects into a specific layer that stored in Oracle database.
Firstly I set a new QgsFeature(), then set geometry and attributes. Finally perform addFeatures() action.

nodes_pnt = QgsProject.instance().mapLayersByName(u'PLOTPOINT')[0]
layer = iface.activeLayer()
out_geom = layer.selectedFeatures()[0].geometry()
feature_geom = QgsFeature()
feature_geom.setFields(nodes_pnt.fields())
feature_geom.setGeometry(out_geom)
feature_geom['DESCRIPTION_'] = "new_object"
nodes_pnt.dataProvider().addFeatures([feature_geom])

But instead of adding new feature to the layer, QGIS throws an error like

SQL: ORA-01400: cannot insert NULL into ("GIS"."PLOTPOINT"."DOC_ID")

I tried another way without setting fields into feature_geom.

feature_geom.setGeometry(out_geom)
feature_geom.setAttributes(['test', "", "", "","","","",""])
nodes_pnt.dataProvider().addFeatures([feature_geom])

However I still got the same error.

I thought about direct connection to the Oracle database using QSqlDatabase to insert values with SQL code, but here I even did not got any connection though parameters were right.

from qgis.PyQt.QtSql import QSqlDatabase, QSqlQuery
db = QSqlDatabase("QOCI")
db.setHostName("test.local:1521\oradb")
db.setPort(1521)
db.setDatabaseName("oradb")
db.setUserName("gis")
db.setPassword("gis")
db.open()

Here I had my db.open False everytime.

What can be done here?


texts: gdalwarp DEM with different resolution extents

texts:
I'm trying to convert DEM source (LIDAR is original source, in ASC format) to a custom internal format  that has specific resolution and alignment values.  I also want to generate different output resolutions (each lower resolution is 1/9th resolution of higher res).
I've written a GDAL plugin to convert from source DEM to the internal format, but it requires the source data to already be in the correct resolution and alignment.  So I'm trying to figure out the right combination of gdalwarp, gdal_translate and gdal_retile to achieve this.  I need the pixel spacing (resolution) to match the target format resolution and I need those pixels to align with the target formats restrictions.  I'm fine cutting data off the edges of the source.
I use python to calculate the correct commands and run them.  I end up with (target SRS is WGS84/EGM96 - hopefully I got the t_srs right).
gdalwarp -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=./egm96_15.gtx" -tap -et 0 -r bilinear -tr <res> <res> -te <res_xmin> <res_ymin> <res_xmax> <res_ymax> <src.vrt> <dest.tif>

where <res> is specific to the resolution (such as 0.00001028806584362140 for 1/27th AS).  The -te values are specific to the alignment requirements for that specific resolution.  So I figure out the coordinates of the target format file that is "inside" the source data (cutting off the edge), and make sure the extent is a multiple of the target file size (again staying inside the bounds of the source data).
Then I use gdal_retile to split up <dest.tif> into individual files to fit what my custom format needs (300x300, 1 pixel overlap, so I do 301x301).
Finally I use gdal_translate with my custom plugin to convert the split files to my custom format.
This seems to work for my different output resolutions, but when I try and view the data, the higher res data does not line up with the lower res data - when I load the higher res data over the top of the higher res data, terrain features do not line up as expected.
So first question: when gdalwarp process -te does it try and fit the source data to that extent (i.e. 'squashing' the source data to fit, assuming -te is smaller than source extent), or does it 'cut off' the part of the source data that is outside that extent?  I'm trying to get it to cut it off.  If it's 'squashing' instead, that might be why I'm seeing alignment problems.
[edit : comment indicates that gdalwarp "cuts" the data]
I've also tried:
gdalwarp -t_srs "..." -tr <res> <res> <src.vrt> <temp.tif>
gdal_translate <res_xoff> <res_yoff> <res_xsize> <res_ysize> <temp.tif> <dest.tif>

i.e. just have gdalwap reproject without trying to "trim" or "align" the output, then use gdal_translate to do the trimming to match my desired alignment.  But this seems to end up with the same problem, and I'm concerned that the pixel alignment won't match up.
Does this sound like I'm going about this the right way?  Any suggestions to try and figure how to re-project and trim the data to match resolution and alignment restrictions?
Thanks!


texts: Extract time series from ImageCollection for multiple points - Google Earth Engine

texts:
I have the Sentinel 2 (COPERNICUS/S2) ImageCollection and a FeatureCollection of points. I want to extract a time series for each of the points from the band information in the ImageCollection.

So far I have the following code:

// Create a polygon to subset sentinel data
var poly = ee.Geometry.Polygon(
        [[[13.932985323828689, -14.515258012418162],
          [13.932985323828689, -16.041324207881946],
          [15.690797823828689, -16.041324207881946],
          [15.690797823828689, -14.515258012418162]]], null, false);

// Import sentinel data
var sentinel2 = ee.ImageCollection("COPERNICUS/S2")
  .filterBounds(poly)  // TESTING - Filter by polygon
  .filterDate('2016-02-01', '2019-02-10');  // TESTING - Filter on date

// Create FeatureCollection of points
var features = [
  ee.Feature(ee.Geometry.Point(14.5, -15)),
  ee.Feature(ee.Geometry.Point(14.8, -15.1)),
  ee.Feature(ee.Geometry.Point(14.9, -15.2))
  ];

var points = ee.FeatureCollection(features);

// Define a function that works on a single point in the FeatureCollection
var pointExtract = function(feat) {

  // Get geometry of point
  var geom = feat.geometry();

  // Define a function that gets a point value for an image
  var pointGet = function(img, f) {
    // Duplicate feature
    var newf = ee.Feature(f);

    // Extract band value for pixel
    var value = img.reduceRegion(ee.Reducer.first(), geom, 30).get('B2');

    // Return as a string in the duplicated feature
    return ee.Feature(
      newf.set(ee.String(value)));
  };

  // Iterate over imageCollection to get time series
  var newfeat = ee.Feature(sentinel2.iterate(pointGet, feat));
  return newfeat;
};

// Iterate over points in FeatureCollection to extract time series
var output = points.map(pointExtract);

// Inspect output
print(output);

But I get the error:

FeatureCollection (Error)
Error in map(ID=2):
String: Parameter 'input' is required.


texts: Splitting polygon into three equal parts using NetTopologySuite and C#?

texts:
I have a polygon, I want to split it into three equal parts, I have a scenario as below:

User must select a vertex on polygon
User must add direction to start split on
Based on the node and direction we must select a line segment on polygon
??????

I have no idea how to go farther and then do what to split that line segment. In fact here are some problems:

How we can find if the line segment is dividable into for example 3 parts?
how we van determine the length of each segments on line to divide?
If line segment changes its direction what will happen?

Can someone give me some hints how to reach a flowchart on this?

I need the methods to do it,algorithm and things like that, I use .netTopologysuit in C# .net


texts: Drawing Google Fusion Tables on Esri Silverlight Map?

texts:
I've been asked to add drawing support for Google Fusion Tables to an Esri Silverlight map application.  For example this one showing Full-Year Stream Temperature Monitoring Sites for the Western US.

Ideally, I'd like to be able to handle this just like any other mapservice layer by providing the url and add it to the map.

Does anyone know of an easy way to do this?

A related question is: would doing this be a violation of the Google Fusion Table Terms of Service?

  5.3. You agree not to access (or attempt to access) any of the
  Services by any means other than through the interface that is
  provided by Google, unless you have been specifically allowed to do so
  in a separate agreement with Google.

Does this mean I have to use Googles javascript API method if I want to show layers from fusion tables?


texts: PDAL works in Windows command but not in OSGeo4WShell?

texts:
I want to import a pointcloud into a PostGIS database. AFAIK you accomplish this with PDAL which comes with OSGeo4W.

I was following this tutorial: OSGeo4W LiDAR tutorial
and also this tutorial PDAL.io tutorial.
Both of them use the osgeo4w Shell, but for me the command
    pdal
throws an error in the OSGeo4W Shell:

  Command is unknown or can't be found.

When I execute any command in cmd.exe the pdal command seems to work.

Can anyone explain why this is?

I want to understand this properly to avoid any mistakes.


texts: How to customize inputs in Iterate Feature Class - ModelBuilder 10.1

texts:
I am using the Iterate Feature Class tool in ModelBuilder to iterate through feature classes within a Geodatabase. I would like to add the capability to custom select which feature classes within a single GDB to iterate through. This desired capability is very similar to the "Input Datasets" step of Merge (Data Managment).

The Iterate Feature Class tool only provides you with 3 ways to filter:

-Feature Type (line, polygon, etc)
-Wildcard
-Recursive

The ability to custom select which feature classes to iterate does not exist, can someone please help me figure out how to customize the input? Thanks!


texts: .dwg CAD file has annotations but they're not showing up in ArcMap

texts:
I have a CAD file with tree data I need to put into ArcMap. I need the species information of each tree point. When I open up the .dwg file in Autodesk TrueView 2018, all the necessary shows up as can be seen in this image:

When I import them into ArcMap as a .dwg file (without converting into a geodatabase format) I don't see the annotations (as can be seen in the following image). Also, these points ended up becoming one huge conjoined polyline rather than individual polylines or points, and there's no species attributes attached to the layer.

So I thought I just had to convert it with the "CAD to geodatabase" tool. However, when I did that what I got was different from the .dwg file, with 98% of my tree data missing except for a few random polygons here and there.

ArcMap did warn me of some errors while the tool was executing. Does this have anything to do with why my annotations are not appearing?

I am able to extract every tree point from that one conjoined polyline. What I'm having trouble with is getting the annotation over each tree point so that I can identify each tree.


texts: How do I add OpenLayers “overlays” using Vue.js

texts:
I can't figure out how to implement an OpenLayers "overlay" in my Vue.js app. There doesn't seem to be much guidance online about how to do this. The Vue component below gives me the map with the vectorlayer, but the overlay does not appear.
PS: Also, if you think there are easier ways of implementing OpenLayers into a Vue framework recommendations are highly welcome. This approach doesn't feel like it's the best.
<template>

</template>


texts: Looping through two raster datasets and performing raster calculation on them with Arcpy

texts:
I have 2 raster datasets of temperature and precipitation from WorldClim. Is there any way to loop through 2 raster datasets at the same time and perform raster calculation on them. So that the combination on every loop is something like this:
temp1 and prec1
temp2 and prec2
..
temp12 and prec 12

I know how to loop through 1 raster dataset and perform raster calculation on it which looks like this:
import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace = r'D:\Data'
rasters = arcpy.ListRasters('*.tif*')

for raster in rasters:
    outraster = Raster(raster) > 0.1
    #Save temp raster to disk with new name
    outraster.save(raster.replace('.tif', '_temp.tif'))

print('Done Processing')


texts:
How to input rpc text file in gdalwarp taransformation. I tought as follow

C:\Users\admin>gdalwarp -t_srs EPSG:4326 -rpc input_rpc.txt  input.tif output.tif

Its giving error as input_rpc.txt not recognized as supported file format


texts: Coordinate background colour using Print composer in QGIS

texts:
How to add coordinate background in Print composer?
Example in attached picture with coordinates, which is not quite visible.


texts: Getting all users in ArcGIS Online organization?

texts:
I want to return a list of all users in my ArcGIS Online (AGOL) organization. However, the result returned in gis.users.search() appears limited to 100 users.

Is there any way to return all users or to perform multiple searches until all users are found?


texts: Python crashing and receiving "EOFError: [Errno 10054]" when encountering arcpy.Exists()

texts:
I have written a script that looks at layer (.lyr) files inside of multiple directories and successfully repoints their corresponding datasets to a new ArcSDE connection. When I test my script on a sub-directory containing a bunch of rasters, upon encountering arcpy.Exists(lyr.datasetName), my script will crash on random "File System Rasters," that is, rasters whose datasets are apparently stored in folders, not in ArcSDE/SQL Server. Since these datasets are not stored in ArcSDE, the arcpy.Exists() function should simply return a value of False and pass over these rasters. And it does successfully pass over most of them, but then it will crash on others and not necessarily the same dataset each time.

In order, I will receive the following messages:

After clicking through the messages, I see the following in my interpreter:

LiDAR Post-Sandy DEM, 1-m grid (USACE) KG Test.lyr
Layer is a group layer

lyr: 2012_PostSandy_NY_NJ_000248_1mGrid.tif
2012_PostSandy_NY_NJ_000248_1mGrid.tif
Traceback (most recent call last):
  File "<string>", line 73, in execInThread
  File "C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\netref.py", line 196, in __call__
  File "C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\netref.py", line 71, in syncreq
  File "C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\protocol.py", line 431, in sync_request
  File "C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\protocol.py", line 379, in serve
  File "C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\protocol.py", line 337, in _recv
  File "C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\channel.py", line 50, in recv
  File "C:\Program Files (x86)\PyScripter\Lib\rpyc.zip\rpyc\core\stream.py", line 166, in read
EOFError: [Errno 10054] An existing connection was forcibly closed by the remote host

Here is a Stack Exchange-friendly version of my code:

import arcpy, os

layersFolder = r'MY_LAYERS_PATH'
arcpy.env.workspace = r'MY_SDE_CONNECTION'
newSDEPath = arcpy.env.workspace

walk = arcpy.da.Walk(layersFolder, datatype="Layer")
for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        print filename
        path = os.path.join(dirpath, filename)
        layer = arcpy.mapping.Layer(path)

        layerList = []

        if layer.isGroupLayer:
            for sublayer in layer:
                layerList.append(sublayer)
        else: layerList.append(layer)

        for singleLayer in layerList:
            lyr = singleLayer
            print "lyr: {}".format(lyr)

            if lyr.supports("dataSource") and lyr.supports("datasetName"):
                print lyr.datasetName

                # This is where Python is crashing
                if arcpy.Exists(lyr.datasetName):

                    if lyr.isFeatureLayer == True:
                        print "Feature layer"

                    elif lyr.isRasterLayer == True:
                        print "Raster layer"

                    else:
                        print "PASSING LAYER: {}\n".format(lyr.name)
                        # break out of current loop
                        break

                    lyr.findAndReplaceWorkspacePath(lyr.workspacePath, newSDEPath, True)

                else:
                    print "PASSING LAYER: {}\n".format(lyr.name)

        layer.save()

print "Script Complete."

UPDATE:
To get around this issue, I've added the following line of code:

if "sde.SDE" in lyr.dataSource:

before I use:

if arcpy.Exists(lyr.datasetName):

I think it's a good temporary workaround, but I'd still like to know why Python is crashing...


texts:
How to skip missing pixel values due they are masked? In the code below, empty lines with missing values in the spectral bands appear in the exported CSV file.

var geometry = ee.Geometry.MultiPoint(
        [[-77.10757333516824,34.93299957943245],
         [-77.05352362680463,35.141503960507634]]);

// Function to mask clouds based on the pixel_qa band of Landsat 8 SR data.
function maskL8sr(image) {
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  var qa = image.select('pixel_qa');
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

var pixelsize = 30;
var data = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
              .filterDate('2016-01-01', '2016-03-01')
              .filterBounds(geometry)
              .map(maskL8sr)
              .getRegion(geometry, pixelsize);

var keys = ee.List(data.get(0));
var featureCollection = ee.FeatureCollection(data.slice(1).map(function(singleData){
  singleData = ee.List(singleData);
  var dict = ee.Dictionary.fromLists(keys, singleData);
  var point = ee.Geometry.Point([dict.get('longitude'), dict.get('latitude')]);
  return ee.Feature(point, dict);
}));

Export.table.toDrive({
  collection: featureCollection,
  description: 'CSV'
});

Link to code


texts: Converting a shapefile from easting/northing to lat/long

texts:
I am working with an open data shapefile that uses easting/northing coordinates.  My geocoded data set uses lat/long, so I am trying to reproject/transform the shapefile to use lat/long.

```{r}
nwk_shp <- st_read('http://data.ci.newark.nj.us/dataset/61038922-a789-4abe-9928-79f9fe96c307/resource/90ead04f-dffd-4b3a-93c7-39cd4b190d4f/download/superneighborhoodswgs84.geojson') %>%
  as('Spatial')
```

Reading layer `OGRGeoJSON' from data source `http://data.ci.newark.nj.us/dataset/61038922-a789-4abe-9928-79f9fe96c307/resource/90ead04f-dffd-4b3a-93c7-39cd4b190d4f/download/superneighborhoodswgs84.geojson' using driver `GeoJSON'
Simple feature collection with 21 features and 8 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 561009.6 ymin: 671183.6 xmax: 598979.6 ymax: 712226.2
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

The shapefile uses northing/easting coordinates:

```{r}
nwk_shp@data[1:3, ] %>% select(-AREA, -PERIMETER, -CATEGORY)
```
1   0     18         LOWER CLINTON HILL   686419.4  575672.5
2   1      8 MT PLEASANT/LOWER BROADWAY   701184.5  583709.9
3   2     12            LOWER ROSEVILLE   701167.2  577849.7

Here is the projection string on my shapefile:

```{r}
proj4string(nwk_shp)
```
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Based on this tutorial, and several others, I have tried several variations on spTransform - none of which seem to transform my data into expected lat/longitude coordinates.

For instance:

```{r}

nwk_shp_transformed <- spTransform(nwk_shp, CRS("+proj=longlat +datum=WGS84"))

nwk_shp_transformed@data[1:3, ] %>% select(-AREA, -PERIMETER, -CATEGORY)
```
 FID MSLINK                       NAME NORTHING_Y EASTING_X
1   0     18         LOWER CLINTON HILL   686419.4  575672.5
2   1      8 MT PLEASANT/LOWER BROADWAY   701184.5  583709.9
3   2     12            LOWER ROSEVILLE   701167.2  577849.7

Depending on what common lat/long CRS I use, I get a 'non-finite transformation detected':

```{r}

nwk_shp_again <- spTransform(nwk_shp, CRS("+init=epsg:4238"))

nwk_shp_again@data[1:3, ] %>% select(-AREA, -PERIMETER, -CATEGORY)
```
35 projected point(s) not finitenon finite transformation detected:
[1] 574389.9 685357.4      Inf      Inf
Error in .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  failure in Polygons 1 Polygon 1 points 1

All of this leads me to suspect that there's something about either

the proj string of the original shapefile?
something omitted from the CRS I am using in my spTransform command?
at work here.

I've read the 'Overview of Coordinate Reference Systems in R' but am still stuck.

Desired output would be a shapefile / SpatialPolygonsDataFrame with coordinates in lat / longitude - eg, Latitudes around 40 degrees, and Longitudes around -74 degrees.


texts: How to invert a Buffer(get the inside of a ring buffer)?

texts:
I have a line layer of a road network around which I create a buffer zone of 1000 m with Analysis Tools>Proximity>Buffer. But I am
interested in the areas outside the buffer zone, and not the buffer zones themselves.

Is there a way to invert the buffer, or which is the best way to pick out the the polygons outside the buffer area?


texts: PostGIS: incorrect interpretation of a polygon that intersects the 180th meridian

texts:
I try to find out all geopoints which intersect the polygon set as a parameter.

The problem is when I pass polygon that roughly covers an area of Bering Strait (nearby 180 longitude):

So I use the query:

SELECT ST_AsText(l.geo_point)
FROM "lightnings" "l"
WHERE (ST_Intersects(ST_GeomFromText('Polygon((132.0 40.0, -148.0 40.0, -148.0 -8.0, 132.0 -8.0, 132.0 40.0) )', 4326), geo_point));

As you can see, vertexes are set in a correct order, clockwise, from the North-West. But the result covers outside area and including whole other world.

For example, in result:

POINT(75.5637 40.0434)

The problem doesn't touch the 0th meridian.

Elementary test:

SELECT ST_Area(ST_GeomFromText('Polygon((0.0 60.0, 10.0 60.0, 10.0 40.0, 0.0 40.0, 0.0 60.0) )', 4326))
UNION ALL
SELECT ST_Area(ST_GeomFromText('Polygon((-5.0 60.0, 5.0 60.0, 5.0 40.0, -5.0 40.0, -5.0 60.0) )', 4326))
UNION ALL
SELECT ST_Area(ST_GeomFromText('Polygon((175.0 60.0, -175.0 60.0, -175.0 40.0, 175.0 40.0, 175.0 60.0) )', 4326))

gives the result:

200
200
7000

Is there any simple trick to force PostGIS understanding me? I don't like an idea to divide the polygon...


texts:
I ran 10 classifications (Random Trees) in ArcMap. Each time the output will be somewhat different.
I have made a final map based on the majority value of each pixel using cell statistics (majority).
I have also produced a map where the number of different class values is reported using cell statistics (variety).

Now I want a map that reports the count of the majority value. E.g. the majority value was 2, how many times out of the 10 did value 2 occur for each pixel.

I know how I could do this using several steps, but as I am running the iteration of 10 classification on 8 different areas this would be extremely time consuming.

Is there a tool or a simple way to estimate or count the number of times the majority value was reported for each cell?


texts: Creating functional spatial index of ST_Geometry Feature Class in ArcSDE for Oracle?

texts: Copied an SDE feature class (with blob(sdebinary(?) geom) into the same SDE GDB (Oracle 11g/SDE 9.3.1) but with the ST_geometry configuration keyword.  I can see that the feature class was created, all the fields come across fine, I can query it with st_geom functions etc., BUT in ArcMap, ArcCatalog, the features don't draw. (*note that when selected I can see the outline of the shape just fine)  When I remove the spatial index the features draw, add the spatial index back, poof no draw.  It would seem that ArcGIS cannot create a functional spatial index of my ST_Geometry features.


texts: Understanding hydrology and DEM results?

texts:
In first picture, you can see an image of a lake(Bellandur lake, bengaluru, india). In second picture, you can see the flow accumulation data for the same area I made(after "fill"ing the sinks). As you can see, the water accumulates rightward of the lake where there are houses and city.

Why is this happening?

Instead of flowing rightward, shouldn't the water flow inward the lake?


texts: How to easily enter a point feature in ArcGIS? (arcmap10)

texts:
My job requires me to enter point features into ArcGIS 10 which correlate to underground RFID markers that mark underground telecom infrastructure. Right now, the workflow requires someone to go out, find a marker, and then draw sketch with measurements to objects on the ground that are visible from an aerial photo to use as reference points. Then, I take measurements off those same reference points in ArcGIS to find the spot where I need to put the point feature. I use the normal measuring tool for this, which has a serious drawback. I've found no way to keep my measurement lines on the screen (or drop a temporary marker or something) when I switch to the edit feature tool to actually drop in the point. Since the points are often out in the middle of streets or fields, I can't just "eyeball" it. This process is too imprecise already.

The workaround I've come up with is using a sticky note and sticking it on my monitor with the corner pointing at the spot where the point should go before I switch tools. There has to be a better way. Any suggestions?

I'd also be interested to hear if there are workflows out there using GPS devices like the Trimble Geoexplorer 6000 series which people have found to be fast enough to justify the investment in the equipment. What we're doing now is labor intensive and error prone. Help me convince my boss to drop $5k on a sub-decimeter GPS!

Thanks!


texts:
The Web AppBuilder is a pure dojo, and Node.js application to build webmap applications by analyst that are not developers. However there is also a developer version where you can extend or build your own custom widgets.

In the Esri example: https://developers.arcgis.com/web-appbuilder/guide/create-a-feature-action-in-your-widget.htm they provide a code snippet but I'm a bit discombabulated. I can see dojo/_base/declare in the define module, but the jimu/BaseFeatureAction,
  jimu/WidgetManager, what js library it is?  What is jimu.js and where can I find more information on it?

define([
  'dojo/_base/declare',
  'jimu/BaseFeatureAction',
  'jimu/WidgetManager'
], function(declare, BaseFeatureAction, WidgetManager){
  var clazz = declare(BaseFeatureAction, {

    iconFormat: 'png',

    isFeatureSupported: function(featureSet){
      return featureSet.features.length > 0 && featureSet.features[0].geometry.type !== 'point';
    },

    onExecute: function(featureSet){
      WidgetManager.getInstance().triggerWidgetOpen(this.widgetId)
      .then(function(myWidget) {
        var vertexCount = 0;
        featureSet.features.forEach(function(f){
          f.geometry.rings.forEach(function(r){
            vertexCount += r.length;
          });
        });
        myWidget.showVertexCount(vertexCount);
      });
    }

  });
  return clazz;
});


texts: Mapinfo 10.5 and 8.5 compatibility? Unable to open TAB file

texts:
I have recently submitted some data to a client as a TAB file. They are using v8.5 and I am using v10.5, the data is just one .TAB (with the associated other files). And all the field types are character, so I thought this wouldn't be an issue. Are there any other things I may have overlooked when saving my file?

Regards,
Alan


texts: In GDAL are there equivalent vector functions for gpkg files similar to those provided by GetRasterBand?

texts:
I am using GDAL (C++) to read raster files into a buffer for example:
cplErr = resDS->GetRasterBand(1)->RasterIO(GF_Read, xPixelUL, yPixelUL, pWidth, pHeight, buff3, pWidth, pHeight, GDT_Float32, 0, 0);

as well as for getting the pixel distribution:
rcutDS->GetRasterBand(1)->GetHistogram(mapData.sCat - 0.5, mapData.eCat + 0.5, mapData.nCat, histogram, FALSE, FALSE, GDALDummyProgress, nullptr);

Additionally, I have vector maps in gpkg format. Converting the vector files to raster files doesn't seem to be a good option. Higher resolution conversions have proven to be difficult.
Are there equivalent methods to the ones mentioned above for dealing with vector files such as gpkg in GDAL? GDAL would be best, but if that's not an option I'm open to exploring other well maintained C++ libraries.


texts: No Shapely geometry can be created from null value during dissolve in Geopandas?

texts:
d1= gpd.read_file(aq.shp)
asttom = gpd.read_file(es.shp)
data2 = d1[d1['TYPE'].str.startswith(('01'))]
dslv=data2.dissolve(by='TYPE')

While it reaches the dissolve procedure it causes the error. How to work with it?


texts: Set the geographic coordinates of a raster map on QGIS using two points

texts:
I am making a map from a raster that obviously has not geographic coordinates. I plan to eventually put my digital map on a world map like OpenStreepMap with the correct scale and positions.

I already know the exact coordinates where this map should appear. I know where my map would fit in.

I was wondering if it was possible to adjust QGIS coordinate system and orientation by using two points for which I would provide the geographic coordinates. Also the map is not oriented towards the north but has an angle of 11° 40' from it, this is why I would want to use two points, so the grid would also have the correct orientation.


texts: How to recategorize features in GRASS by geometry, like the ArcGIS Sort tool?

texts: The ArcGIS Sort tool creates a new feature class with features sorted by geometry. For example, it can number features IDs from north to south, etc. Is there a way to renumber the cat field in a GRASS vector map based on geometry in a similar way?


texts: Hillshade input field for z factor keeps rounding value in QGIS?

texts:
In QGIS 3.8.0, the Hillshade tool
when I try to input the value for the z-factor which for my case should be
0.00000956 exactly, after you remove focus of the input field QGIS rounds that number to 0.00001

This goes for other tools as well not just Hillshade and it was working just fine as expected before, so I guess this some kind of a configuration issue maybe ? I uninstalled and reinstalled but it is still there.


texts: Creating Polygon Shapefile from Point data using QGIS?

texts:
I have a point shapefile that has data for application rate, distance travelled, and swath width.  In another software package it can show the points as shown.

I know how to make QGIS show the point data this way.  However, I want to be able to use the distance travelled and swath information to create a polygon surrounding each point and assign the rate to that polygon.  In the end I want a map that looks like this:

My other software package will display this map, but there is no way for me to export it as a shapefile so I can use it in another program.  Does anyone know of a way I can get QGIS to create this type of map?  I've searched around but couldn't find anything that matches what I'm trying to do.

I've worked with the software company that created the map and they do not have a way currently to export that layer.  I'll see if I can get them to add that in the future, however it of course doesn't help me right now.

I've loaded a sample file with the data and information I want to use at:

Sample_Data

The data in the table I want to use is:  "Swth_Wdth_" and "Distance_f" (both values are in feet) which create the rectangle around the point.  I want to assign the rest of the attributes for each point to the polygon that is created (however "Rt_Apd_Liq" is the most important).  If it helps, "Track_deg_" is the heading value (direction) when the point was recorded.


texts: Is there a way to add system date or time in print composer?

texts: Is there a way in the Print Composer of QGIS to create a label that will automatically update the date (systems date) instead of having a label that you have to manually have to update?


texts: I digitized tracks to the scale 1/4000 and I want to get out a map at 1/50000. I need a tool in ArcGIS or QGIS with its parameters that will help me to generalize the tracks to the scale 1/50000


texts:
I am new to GIS.  I downloaded GRASS and the sample data sets and am attempting to use the tutorial at: http://grasswiki.osgeo.org/wiki/Quick_wxGUI_tutorial.  During the step, "Create a random surface", I get the error:

ERROR: G_malloc: unable to allocate 536870912 bytes of memory at spec_syn.c:50

My computer should be plenty powerful...
8GB RAM, Windows 8, 64-Bit

Any ideas why I can't seem to get through even the most basic tutorial?


texts: PostGIS - "parse error - invalid geometry" after using ST_Multi - but ST_IsValid returns true

texts:
I have a database with OpenStreetMap polygons, points and relations in separate tables. Now I wanted to build selection queries for the different features and created schemas for each zoom level (5k, 10k, 20k and so on). However, I'm unable to insert multipolygons into the "landuse_background" table, even though the geometry column is configured for multipolygons:

CREATE SCHEMA "5k";
CREATE EXTENSION IF NOT EXISTS postgis WITH SCHEMA "5k";
CREATE TABLE "5k"."landuse_background" (
    "id" BIGINT NOT NULL UNIQUE PRIMARY KEY,
    "geom" geometry(MULTIPOLYGON, 4326) NOT NULL,
    "name" TEXT
);

As a test, I wanted to insert every polygon as a multipolygon (so that I can later on combine them):

INSERT INTO "5k"."landuse_background" (
    SELECT osm_id as id, name, ST_Multi(geom)
    FROM "public"."osm_polygon"
    GROUP BY osm_id
);

However, I get the error:

ERROR: parse error - invalid geometry

If I only execute the "select" part of my query, I get valid polygons:

SELECT
    osm_id as id,
    name,
    ST_AsText(ST_Multi(geom)) as geom,
    ST_IsValidReason(ST_Multi(geom)) as is_valid
FROM public.osm_polygon
GROUP BY osm_id;

I get valid geometry:

8023823         "Schlossgarten"             "MULTIPOLYGON(((8.3999983 49.0136811,8.3999971 49.0136871,8.4000414 49.0136908, (...))))"   "Valid Geometry"
8085105         "Alter Flugplatz"           "MULTIPOLYGON(((8.3824917 49.0242842,8.3820936 49.0217043,8.3821084 49.0216129, (...))))"   "Valid Geometry"
19782878        "Hermann-Ehlers-Kolleg"     "MULTIPOLYGON(((8.3971125 49.0156717,8.3971956 49.015771,8.3971811 49.0157767, (...))))"    "Valid Geometry"
20909486        "Augustin-Bea-Haus"         "MULTIPOLYGON(((8.3958855 49.0162626,8.3958992 49.0162606,8.3958896 49.0162402, (...))))"   "Valid Geometry"
20911412        "Reinhold-Schneider-Haus"   "MULTIPOLYGON(((8.3965158 49.0161899,8.396547 49.016189,8.3965608 49.0162083, (...))))"     "Valid Geometry"
23138236        "Domus 7"                   "MULTIPOLYGON(((8.3952316 49.0161717,8.3952569 49.016172,8.3952871 49.0161686, (...))))"    "Valid Geometry"

ALL geometries are valid, selecting invalid polygons returns zero rows.

So I thought, it might be an issue with the CREATE TABLE statement and checked the column types of the landuse_background table:

SELECT * FROM information_schema.columns
WHERE table_schema = '5k'
AND table_name   = 'landuse_background';

table_schema    table_name              column_name     data_type       udt_schema     udt_name

"5k"            "landuse_background"    "id"            "bigint"        "pg_catalog"   "int8"
"5k"            "landuse_background"    "geom"          "USER-DEFINED"  "public"       "geometry"
"5k"            "landuse_background"    "name"          "text"          "pg_catalog"   "text"

The only noticable thing is that the geom column has a data type of USER-DEFINED instead of the expected MultiPolygon. I've tried several identifiers in the CREATE TABLE statement, such as MultiPolygon, MULTIPOLYGON, "MultiPolygon", but none seem to work.

I've already looked at the similar answers

PostGIS Update parse error – invalid geometry
ERROR: parse error - invalid geometry (Postgis)
parse error - invalid geometry even though ST_ISVALID returns true

... but they don't match my problem.

PostgreSQL version: 9.6.6, PostGIS version: 2.3.1 r15264


texts:
I have a polygon which I want to split using overlapping polygons in order to create 3 polygons (with the original attributes)
How can i achieve it using R ?


texts:
I am using Ordnance Survey's OpenMap - Local in Mapnik 3.0.
Ordnance Survey make stylesheets available here:
https://github.com/OrdnanceSurvey/OS-OpenMap-Local-stylesheets
Is there a way to convert any of these formats (ESRI, GML or Geopackage) into the XML I need for Mapnik?
Or alternatively, is there a way to get Mapnik to read any of these formats directly?


texts: How to correctly merge the value in the columns, in the case of using the widget type the Value map, with QGIS

texts:
Trying to do work that was previously done in ArcGIS. There is a table when filling 4 columns, these values are added in the Sum column.

In QGIS, I decided to use the Value Map for fields

For the Sum field, I use the following expression

But as a result, I lose zeros
before number

Where did I go wrong and how do I get values with zeros when field type  is real and SUM field is string?


texts:
I am exporting contour lines from a raster file which was obtained by clipping a shapefile to a much larger raster file. The contour lines obtained are square-ish as in the second image below.
My first attempt in generating contour lines looks like this. The lines are smooth enough.
And my second attempt looks like the image below. Lines are squared as if they are following an invisible grid.
First time I didn't really pay attention to what I was doing and followed a tutorial or two.
Second time I discovered that projection matters and fiddled with it :). Second screenshot is in ESPG:27700 projection (if this matters, I don't know).
Did you see this before and if yes, how do I fix this? If not, can you point me please in the right direction? What to read, or to what this issue might be related to.
My guess is that when I am clipping, resolution is lost. But then how did I manage to obtain smooth lines first time?


texts: How to find a single point furthest from any road?

texts:
I have a road layer of a country. What method or plugin would you suggest, to find a single point furthest from any road

(but as you can see, there is a coastline, so points in the sea can't count)

The road layer looks like this:


texts: How can we change the symbol to be in right side and the label description will be in the left side in qgis 3.6?


texts: Filtering number of features in GeoServer using SLD?

texts:
I have published point file in geoserver  from postgis table as store. I had huge data in my table(rows) and the data will keep udate daily (at least one point for one day). Now I want to show only the last 10 points added.

How do I restrict the data using Geoserver SLD Function?


texts: Where is the OSGeo4W setup for QGIS 2.6.1?

texts:
The downloads site for QGIS appears to only contain standalone installers. I want to install QGIS 2.6.1 using the OSGeo4W Network Installer (64 bit). When I do this currently, I install QGIS 2.7 which is not what I want.

Can I do this or must I download the OSGeo4W installer for QGIS 2.7 and then install the standalone for 2.6.1 and make that the default?


texts: Why is Overpass-turbo map empty if there is data?

texts:
I started with overpass and am using the overpass-turbo.eu site. How comes that the map is empty when there is data? In particular, the park marked in blue shows up in the xml file (with the tag leisure="park") but its not marked on the map.

I would like to query for parks in the close range of the user.

SOLUTION:

way[leisure="park"](around:4500,47.0546694,15.4570067); out geom;


texts: Delete ONLY selected features in ArcMap / ArcPy

texts:
I have a layer where features are selected by location in relation to one or more other layers (selection_type is ADD_SELECTION) in a loop.
After that, zero, one, some or all features are selected.

I want to delete only the selected features - if there are some.

But the "Delete Features"-Tool deletes either the selected features, or - if none are selected - all of them.
How can I check, if there are any features selected?
The "Get Count"-Tool works just like the "Delete Features"-Tool,

E.g. Layer X has a total of 100 features:

No features are selected, "Get Count" returns 100
All features are selected, "Get Count" returns 100

I would like to delete only the selected features. If there are none, I don't want to delete anything.


texts: Calculating aspect in flat areas in ASTER DEM with ArcGIS for Desktop?

texts: I have an analyzed aspect surface from ASTER DEM 30m resolution in ArcGIS 10.2. But it shows many aspects in flat area.Normaly it must calculate this are as flat (-1). How can I resolve this problem?


texts: True non-projected 3D grid world modeling

texts:
I did a search on this but I did not find anything much relevant yet. I was wondering how come there is no GIS software that would take a sufficiently accurate 3D surface model of the world (such as based on SRTM's) and put geographic features on it instead of a simplified spheroid.

This would obviate the need for flattening the world using projections, hence remove all distortions, enabling us to make much more accurate measurements. We then wouldn't need to use different projections for different purposes or for different regions of the world.

I'd imagine when we look at the interface of this software, we'd be looking at a curved map, flattened only at the screen. When we zoomed out enough, we could see the entire world as a globe (as an option). The major problem with this approach would be the relatively large size of the 3D surface grid, demanding computation power and a need for developing new libraries for computations.

I think the most of these issues could be solved by getting only the portion of the grid necessary for one's mapping needs.


texts: Is there a programmatic way to remove a user from a license

texts:
I'm new to the Enterprise level of the GIS world and working at this level I find that there are many instances where another user has taken an ArcEditor or ArcInfo license when it's needed for a more important situation that has arisen.

My question has two parts.

Is there a way to remove a user from a license programmatically?
Is there a way to transfer a user from one license-type to another? (i.e., ArcEdit to ArcView)

Programmatically would include Python, VBA, VB, or other.

I am familiar with "Where Have All The Licenses Gone" but it's very inconvenient as sometimes we have users that leave their computers on and then leave for the day. Rather than unplugging their systems and potentially losing their work it would be a lot more efficient to just strip a specific license from them or transfer their current license to another license.


texts: Pansharpening radiance/reflectance of Landsat 8 in GRASS 7?

texts:
first at all, I have to say that I'm new GRASS user.

I would like to make a pansharp in GRASS 7 with reflectance Landsat 8 bands. I have revised some related posts and blogs, and I didn't find out if it is possible to make a pansharp after in reflectance bands.

It seems that i.pansharpen only works with 8-bit data, however, Landsat 8 is delivered in 16-bit.

Is it necessary to do a r.rescale with reflectance Landsat  8 data before pansharpening?


texts: Increasing QGIS 3.2 (Bonna) Render Size on my Macbook Air

texts:
The render size for the project appears to be small in the QGIS 3.2 for MacOS.

Any tips how to increase this/make it wider?


texts: Landsat7 - converting dns to at sensor radiance

texts:
I am trying to convert DNs of Landsat 7 ETM+ bands to at sensor radiance values. I know the formula that I can use in the raster calculator, but I am confused about the values in the Metadata. I am looking for Lmax and Lmin values for all bands.

Are these the same as radiance maximum and minimum values?

Qcal min and max values are there in the metadata but what about the Qcal (not min or max) value from which we need to subtact Qcal min in the algorithm.


texts: How to avoid creating corrupt Shapefiles during Editing?

texts:
I have one of my GIS technicians digitizing some lines in QGIS in shapefile format. I don't know how he did it (and neither does he), but somehow the shapefile became corrupt. It was creating random lines or some of the lines he created would just disappear. I went in into ArcCatalogue to see how it looked in ArcGIS and this is what I saw:

Notice the question mark icon where I should see a shapefile 'line' icon. Obviously ArcCatalogue cannot read this file. Also, a second dbf file seemed to have been created with the '_packed' attached to the end. When I look at the shapefile using Windows explorer, I see that there already is a .dbf for the shapefile 'M3_PRE_SMU_lines_10Apr13_SMC.dbf', so I don't know where this _packed shapefile came from and I can't seem to find anything online that speaks to it.
I tried to add this file into ArcMap and I receieved the following error:

The error is pretty self explanatory...number of shapes does not match the number of records. I just don't know why that is occuring. There doesn't seem to be anything online that explains how this is occurring in QGIS, but I do see a couple repair tools. I actually repaired this myself by just opening QGIS, adding the layer, and then right clicking on the layer and 'saving as' another shapefile. So, I've figured a work around, but I'm hoping to find a solution that will stop this from occurring in the first place.
Thanks,
Mike


texts: Landsat 8 OLI Tasseled Cap Transformation Calculation ArcGIS 10.3

texts:
I need to calculate the TCT greenness for Landsat 8 OLI imagery and having a hard time finding a tool.

All the python scripts I have found do not work and I am not experienced with Python at all.

Can anyone please point me to a step by step tutorial or just personal experience with this?


texts: Average Nearest Neighbor does not have all the options

texts:
I have a point feature class that I used the Average Nearest Neighbor from the ArcGIS software to calculate the pattern.

ArcGIS produces output chart that says whether the data is dispersed, random or clustered. But what if it's uniformed? what if it's regular?

Why are these patterns not in the ArcGIS result and what can I do to tell if my point feature class is uniformed or regular using ArcGIS and the "Average Nearest Neighbor tool"?


texts:
I am an avid user of OsmAnd. Some of the maps I have are now becoming a bit outdated - shops that changed name, paths that no longer exist, new paths that are not yet marked.

How can I update a map that I already have on my system?


texts: How the year affects Area Solar Radiation?

texts:
How choosing a different year in "Area Solar Radiation" affects the calculation of the radiation?
I understand that months change some parameters in the calculation formulas, but what about the year? What's the difference between choosing 2012 or 2013, or some other year?


texts: Where to find or how to create from scratch a shapefile of an obscure city?

texts: I am trying to find the shapefiles of a couple of cities in Mongolia, for example: Sukhbaatar (in Selenge, Mongolia).  I haven't been able to find a source that provides these shapefiles, so I thought maybe I could get the city shapefile from OSM.  Currently I try to export the osm file after typing in Sukhbaatar, Selenge, Mongolia and use osmtogeojson (https://tyrasd.github.io/osmtogeojson/) to get it into geojson form first. However, the features are an empty array.  Since that doesn't work, is there an alternative to obtain the shapefile or geojson of an obscure city?


texts: How to correctly use the filter with empty values in PyQGIS?

texts:
In my script, at the moment of having the layers loaded, I filter all the ToC layers with a unique value for all the layers.

I use the setSubsetString('"FIELD" = \' VALUE \ '') method, but when it does not find any field with that value, it throws me the following error: wkb outbounds.

How can I solve this? I am using QGIS 2.14.


texts: Drawing perpendicular line between point and line layer using ArcMap?

texts: I have downspout points that are currently emptying at the foundation. We need to redirect those into the roadway to be collected by the storm sewer. However we need to check that there is enough vertical elevation change between the base of the house and the curb to slope the pipe appropriately (e.g. if the downspout is 50ft away and the elevation at the foundation is 700ft and at the curb 699ft and I need 18" of cover then this will not work). So I need to construct a line from every single downspout point that connects to the edge of pavement line and get the elevation at each end. This line also needs to meet the edge of pavement line at a 90-degree angle. I have over 4000 downspouts so I'm looking for a way to automate the generation of the line segment between the two in such a way that when I calculate the start and end elevation using my DEM I will know my grade change.


texts: Esri File Geodatabase The parameter is (or has) the wrong kind of geometry

texts:
One of my tools was having problems, during debug I encountered a new error based on FDO error codes:
Case -2147220982
    fKnownErrNum = "The object has been deleted and is no longer valid."

Attempting to edit the feature class in ArcMap:

Which suggests that the data editor (somehow) has inserted a geometry that isn't the expected PolyLineZ into the table.
I tried to remove the bad shape with a tool I wrote but on the line pFeature = pFeatCur.NextFeature I get this error message:

Theoretically it should be impossible to insert a row with an incorrect geometry into an Esri file geodatabase, which has me scratching my head and wondering, but academics aside the data represents capture time and I would like to salvage if possible what I can. So far I have tried Repair Geometry, my own tools (with obvious failure) and the usually reliable Export to XML workspace document - which returns a similar error to editing in ArcMap.
Has anyone experienced a similar issue to this and been able to redeem (even partially) the data?


texts: Qgis 2.8 filtering the ten most populated places for a specified altitude

texts:
I try to filter data combining settlement and altitudes.
I try to filter the 10 most populated places for a specified altitude.

I tried this expression without success:

"altitude">1000 AND ORDER BY POPULATION DESC LIMIT 10


texts: Attribute table empty on a layer filled with data

texts:
I have inserted some data from an .xls file and marked it on a map in ArcMap on a layer named "Kafichi"

However when I right click on the layer and click on "Open Attribute Table" it gives me an empty table even though I have some data filled in.

Can anyone help?

I've tried exporting the data from the "Database$" however that did not work.

I can also assume that this is the same reason that when I do a Point Distance on those Coffee values with another layer it gives me an empty table.

Correct me if I'm wrong.


texts: Where can I get US snowfall data with day resolution?

texts:
I need to visualize the snowfall in the USA for a given time period and city.  So you could pick a given day and get a heatmap showing the amount of snowfall at each place in the USA on that day.  Monthly or seasonal snowfall isn't granular enough for this application.

Does anyone know where I could get data like this?  (Ideally it would already be hooked up to a map so I can explore the data before building the whole thing myself.)


texts: Importing LiDAR .txt in QGIS for map conversion


texts:
I'm working with ArcGIS 10.2.
I'm using the next Python code to show some records of an attribute table:

import arcpy
arcpy.MakeFeatureLayer_management(r"F:ArcGIS\Trab\Model\Lin_150_tr.shp", "lyr5")
arcpy.SelectLayerByAttribute_management("lyr5", "NEW_SELECTION", """FID IN (2,10,200,323,349,384,433)""")

But I'm not able to run my Python script again because appears an error like "lyr5 already exists". So I use a funtion:
arcpy.Delete_management(lyr5)
but it gives me an error NameError: name 'lyr5' is not defined.

What can I do to solve this problem, updating my layers, without loss
my original shapefile called Lin_150_tr.shp?
May I do a copy of this shapefile before to do
SelectLayerByAttribute?


texts: Select by attributes from a series of feature classes

texts:
I have 16 feature classes. What I want to do is to select by attributes each feature class using this statement:

    "Shape_Area" < 1100

Is there a way to do this by batch using arcpy or arcpy.mapping? I can't find the right module for "Select by Attributes".


texts: Adding lat/long coordinates to clicked feature on ArcGIS Online?

texts:
When I click a feature in my ArcGIS Online, I would like to customize the popup to display the latitude and longitude coordinates of the feature along with the fields for the feature that was clicked.

Can this be accomplished with ArcGIS Online?


texts: Is there any way to read raster path in POSTGRES and bring to QGIS screen

texts:
.qlr or .vrt files contain a list of raster files path such as : 'C://raster/NYK.tiff'....'C://raster/CAL.tiff'.... etc...

QGIS displays them by reading the embedded paths.

I would like to store these path in a Postgres table. Is there any way for QGIS to read the paths from such table and to display the raster? Let's note that I don't want to store the raster themselves in the DB, only their path.


texts: Extrapolating with IDW (Feature Classes) in ArcGIS Pro?

texts:
I have been struggling to extrapolate beyond the extent of the point features input into the IDW tool of ArcGIS Pro.

I believe IDW interpolation is the method for me, i just need the process to extend beyond the input points slightly to fit a bigger area.

Specifying a mask in the environment settings seems only to clip the same extent (that is that of the input points), to the mask.


texts: Compute parallel lines along a center line in PostGIS

texts:
I have street (road center line) and building polygon tables in my PostgreSQL database. The sample scenario is as follows:

Problem:

I need to compute parallel lines along street at the intersection of 50 meters buffer around street and nearest building polygon on both sides. The desired output scenario is:

What I have tried:

My approach was:

 1) Generate 50m buffer around street layer
 2) get the intersection of buffer and polygons
 3) Compute the distance
 4) Draw offset curves (parallel lines) at both sides of street layer
 5) Merge both curves to get parallel lines at the intersection

Here is my attempt:

    WITH street_buffer AS (
     SELECT
      street.gid street_id,
      street.geom street_geom,
      ST_Buffer(street.geom, 50, 'endcap=square join=round') geom1,
      building.geom  geom2
     FROM street
     LEFT JOIN building on ST_DWithin(building.geom, street.geom, 50)
     ORDER BY street_id
    ),
    selected_buildings AS (
     SELECT
      street_id,
      street_geom,
      ST_Intersection(geom1, geom2) geom
     FROM street_buffer
    ),
    distance AS (
     SELECT
      street_id,
      street_geom,
      ST_Distance(street_geom, geom) as dist
     FROM selected_buildings
    ),
    curves AS (
     SELECT
      street_id,
      ST_OffsetCurve(ST_LineMerge(street_geom), dist) as curve1,
      ST_OffsetCurve(ST_LineMerge(street_geom), -dist) as curve2
     FROM distance
     ORDER BY street_id
    )
    SELECT
     street_id,
     ST_Union(curve1, curve2) geom
    FROM curves
    ORDER BY street_id

The problem with above code is that it doesn't return parallel lines according to the desired output i.e., parallel lines at all polygons intersections are being generated instead of at the intersection of nearest polygons.

EDIT_1:

The actual output of above code is:

While, in the above output only yellow parallel lines (offset curves to nearest polygons on both sides of street) are required:

Can anyone suggest me how to get the desired output?


texts: Automating ArcGIS share service with Python?

texts:
Right now i am sharing a map service in ArcGIS by following these steps.

File --> Share as --> Service --> then click on overwrite an existing service --> select appropriate map --> Click on Analyse  --> Publish service.

Is it possible to do all this task with a small python code?


texts: How to make zonal statistics model more efficient?

texts: I've made this model, to assess some 500-600 rasters for one area, and spit out one huge table. So far, when testing it with 20-30 rasters as inputs it works, but now that I put in the full amount its taking hours. It appears that it's getting hung up on the merge tool. I thought that the collect values tool would hold the individual tables, and then as a last step merge them together, but it seems like that is not the case (?).  Does anyone have any ideas about how to make this more efficient?


texts: Ways to step through features without Browse Features from data reviewer extension or data driven pages?

texts:
I'm looking for ways to browse through a bunch of points that I have so that I can assign each of them a value. Ideally I would like to be able to use a button on the keyboard to step through each individual point (while keeping the same zoom level) so that I can keep my cursor over my attribute table in order to assign a value to the point. I have used Browse Features before and I know that is an obvious solution but I currently do not have the extension. Data driven pages does work, but it doesn't select the feature so it is easy to get lost where I am when I am in the table when I go to assign a point. I also have to move my cursor to click next page which is slow and clunky and I can't seem to assign a keyboard shortcut to it as well. It also doesn't seem to update the value in the drop down menu to what I have assigned to the name field so I can't keep track of where I am in the table with that either.

Does anyone have any ideas? Being able to stay zoomed in is also important since I need to be able to see the color of the pixel from my imagery that my point is sitting on top of.