texts:
I'm preparing a raster layer of combined resistances of landcover, roads and rail and a second raster layer of 'focal nodes' for a Circuitscape Analysis. I prepped the files in R and they look fine in ArcGIS.

But, when I convert my .tif files to .asc, I get this error and then subsequently, Circuitscape does not run when I input them in the software.

gdal_translate('Prot_Areas_all1.tif', 'Prot_Areas_CS.asc', of='AAIGrid',
FORCE_CELLSIZE=TRUE)

Warning 1: Producing a Golden Surfer style file with DX and DY instead
of CELLSIZE since the input pixels are non-square.  Use the
FORCE_CELLSIZE=TRUE creation option to force use of DX for
even though this will be distorted.  Most ASCII Grid readers
(ArcGIS included) do not support the DX and DY parameters.

Circuitscape throws up an error saying 'Raster headers do not match'.

My R skills are pretty rough, so I don't really understand the exact problem I'm facing either.

This is a part of the code I was using for one of the 'Focal Nodes' Raster file:

 #Converting Protected Areas shape file to .asc for circuitscape
 #First convert to raster and then to .asc
 Prot_Areas <- readOGR(*shapefile from my drive*)
 b<-bbox(Boundary_transform) #retrieving spatial boundaries from our project
 boundary
 r<- raster(nrow = 500, ncol = 500, xmn= b[1,1], xmx= b[1,2], ymn= b[2,1],
 ymx= b[2,2])
 Prot_Areas_transform <- spTransform(Prot_Areas, ll) #ll here is the UTM Zone
 44 projection
 Prot_Areas_raster <- rasterize(Prot_Areas_transform,r)
 plot(Prot_Areas_raster)
 Prot_Areas_raster

 #Recalssify all to one value= 1 (these are the focal nodes- which must be
 equal in value)
 m <- c(1, 15 , 1)
 rclmat <- matrix(m, ncol=3, byrow=TRUE)
 Prot_Areas_rc <- reclassify(Prot_Areas_raster, rclmat)
 plot(Prot_Areas_rc)

 # resample to make the same resolution as Land Use Land cover (this is the UTM Zone 44 proj, and 100X100 resolution
 Prot_Areas_reasmpled<-resample(Prot_Areas_rc, LULC_rc, method="ngb")
 plot(Prot_Areas_reasmpled)
 writeRaster(Prot_Areas_reasmpled, 'Prot_Areas_all1.tif', overwrite=T)

 #Creating the .asc file for Circuitscape
 gdal_translate('Prot_Areas_all1.tif', 'Prot_Areas_CS.asc', of='AAIGrid',
 FORCE_CELLSIZE=TRUE)


texts: ArcGIS Field Editor if/else syntax with python

texts:
I am just starting with Python and ArcGIs so here is my question:

I have the column "DB_Low" (Type "Long") and want to fill a new column by changing values. The new column will be a Float-Type field.

So here is what I tried so far:

if !DB_Low! == 0:

  return 0

else:

  return  !DB_Low! + 2.25

But I get the Error: 000539 SyntaxError, invalid syntax (, line 1)

Can someone tell me the mistake in my code?


texts: How to export map from raster with coordinates to image/pdf file?

texts:
I am using QGIS and I'd like to export a piece of map, scale 1:10000. The area I want to export is wider than the QGIS window, so I can't just use the option "export everything you see on screen", because afterwards I'd need to paste 10 or 12 individual files and bring them together, which is kind of time-consuming.

Is there a  way to export a map using polygon (rectangle in my case) coordinates and to convert it to PDF or some image format?


texts: QGIS 3.8.3 Python Ubuntu 18.04 Error : Could not load source layer for INPUT

texts:
I am running a .py file from terminal in Ubuntu 18.04.

python3 script.py

I am running the given below algorithm of QGIS Desktop 3.8.3 in script.py.

processing.run("gdal:cliprasterbymasklayer", {'INPUT':'/home/sreeraj/Desktop/Program/IMAGE.tif','MASK':'/home/sreeraj/Desktop/Program/VECTOR.geojson','SOURCE_CRS':QgsCoordinateReferenceSystem('EPSG:32643'),'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:32643'),'NODATA':None,'ALPHA_BAND':False,'CROP_TO_CUTLINE':True,'KEEP_RESOLUTION':False,'SET_RESOLUTION':False,'X_RESOLUTION':None,'Y_RESOLUTION':None,'MULTITHREADING':False,'OPTIONS':'','DATA_TYPE':0,'OUTPUT':'/home/sreeraj/Desktop/Program/OUTPUT.tif'})

I am getting the following errror :

Could not load source layer for INPUT: /home/sreeraj/Desktop/Program/IMAGE.tif not found

At the same time, when I open QGIS 3.8.3 in Ubuntu 18.04 and execute this same code in Plugins -> Python Console, it is working perfectly without any errors and is creating the output.

I am wondering what is happening here !!! It is the same code I am executing inside .py file and in Python Console of QGIS, but in the case of .py file, I am getting this error.

This error means that the algorithm is not able to find this .tif file. But, then how this same code is working in QGIS Python Console ? How can I solve this error ?

The exact code in script.py is as given below.

import sys
from osgeo import ogr, gdal
from gdalconst import *
from qgis.core import *
import qgis.utils
from qgis.gui import *
from PyQt5.QtGui import *
from PyQt5.QtCore import Qt
from PyQt5.QtCore import QFileInfo
from PyQt5.QtWidgets import QApplication
app = QApplication(sys.argv)
QgsApplication.setPrefixPath("/usr/share/qgis", True)
qgs = QgsApplication([], False)
sys.path.append('/usr/share/qgis/python/plugins')
from qgis.analysis import QgsNativeAlgorithms
import processing
from processing.core.Processing import Processing
qgs.initQgis()
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
from processing.tools import *

processing.run("gdal:cliprasterbymasklayer", {'INPUT':'/home/sreeraj/Desktop/Program/IMAGE.tif','MASK':'/home/sreeraj/Desktop/Program/VECTOR.geojson','SOURCE_CRS':QgsCoordinateReferenceSystem('EPSG:32643'),'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:32643'),'NODATA':None,'ALPHA_BAND':False,'CROP_TO_CUTLINE':True,'KEEP_RESOLUTION':False,'SET_RESOLUTION':False,'X_RESOLUTION':None,'Y_RESOLUTION':None,'MULTITHREADING':False,'OPTIONS':'','DATA_TYPE':0,'OUTPUT':'/home/sreeraj/Desktop/Program/OUTPUT.tif'})

qgs.exitQgis()
print('Done')

PLEASE NOTE THAT THIS IS NOT AN ISSUE WITH THIS ALGORITHM. WHEN I AM TRYING WITH OTHER ALGORITHMS ALSO, THIS SAME ISSUE IS COMING UP.

I doubt that there is a problem with the path convention starting with '/home/' when using these QGIS Algorithms in TERMINAL. Otherwise, it may be a issue with the initialization of processing and/or importing. So, please check in detail the headers used in this .py file.


texts: I've read a number of articles on RTK GPS.  I'm only familiar with using RTK GPS with an RTK Base Station, and I understand how the corrections can be broadcast to a rover.  When would you use RTK GPS without using a base station, and is this practical?


texts:
I'm trying to import a CSV file in QGIS 3.4 using the common method of "Add Layer" > "Add Delimited Text Layer". But QGIS is reading every row of the dataset as a column title, and nothing is in a row below those columns. See the screenshot below:

I thought there might be some weird formatting issues, so I tried saving the dataset various ways (xlsx, csv, txt) but no success.

Then I tried to create a very simple csv file from scratch as a test, which looks like this:

And when trying to "Add Delimited Text Layer" with this file I get the same problem. All rows are read as column headers:

What obvious thing am I missing? It is as if this tool suddenly does not work, though I have used it many times before.


texts: Arcpy Creating Multi-page Table/PDF with Reportlab/PyPDF2

texts:
I am writing a scrip to batch automate hundreds of map PDF's. Along with the PDF's I generate a unique Table/reportlist for each map from the 'receptors' that fall within each boundary of the map. The amount of receptors and consequently the size of the Table/reportlist can vary from 0 receptors to hundreds. Currently I am generating the Table/receptorlist from a merged clipped receptors for each map and then using the arcpy table to table conversion

PointL = arcpy.TableToTable_conversion(listlayer, gdb, pointsellayer)

I then use 2 arcpy.searchcursors(explain why I use 2 later) to open up and read that table. I then paste the info I need in the order and structure I want onto a ReportLab Canvas.I then use PyPDF2 to merge the map and the receporlist together into one PDF.

My Problem: sometimes as I have previously stated there are hundreds of receptors to copy and write from the Table to the receptor list. What I have been able to do (seems clunky) is write to a one page then use the Show.Page() function to move to the next page using an if statement with a counter.

*skelton outline of my code below(left out all unimportant info for the question)

pointtable = gp.Searchcursor(PointL)
pointtable2 = gp.Searchcursor(PointL)
pdffn = "ReceptorList.pdf"
c = canvas.Canvas(pdffn, pagesize=landscape(legal))
c.drawString(10, 580, "ReceptorID")
loc = 565
z=0
nlist = []
for info in pointtable:
    z+=1     #counter
    recid = info.ID  #example rows
    c.drawString(30, loc, recid)
    c.drawString(other rows...)
    loc = loc-10
    nlist.append(recid)  #storing id in list
    if loc <= 10:
        c.showPage()                    #showPage() = go to next page
        print "ReceptorListpage1"
        loc = 10
        break
if loc == 10:     #page 2 if needed
    loc = 580
    for info2 in pointtable2:
        recid2 = info2.ID
        if recid2 not in nlist:
            z+=1
            c.drawString(other rows...)
            nlist.append(recid2)
    c.showPage()
c.save()
output = PdfFileWriter()
file1 = PdfFileReader(file("outputmap.pdf", "rb")) #created map in this loop
file2 = PdfFileReader(file("ReceptorList.pdf", "rb")) #created receptor pdf         in loop
if z > 64:
    output.addPage(file1.getPage(0)) #map
    output.addPage(file2.getPage(0)) #page 1 list
    output.addPage(file2.getPage(1)) #page 2 list
    outputstream = file("FacilityName+.pdf", "wb") #naming new pdf fac_name
    output.write(outputstream) #merging both pdf's into one!
    outputstream.close()
else:
    output.addPage(file1.getPage(0)) #map
    output.addPage(file2.getPage(0)) #page 1 list
    outputstream = file("FacilityName+.pdf", "wb") #naming new pdf fac_name
    output.write(outputstream) #merging both pdf's into one!
    outputstream.close()

I am able to generate a multi page pdf for a receptor list but there are times when 2 pages is simply not enough and if i were to write more code to for more PDF pages for a receptor it would be very difficult to write because of all my counters, IF statements. Is there another way I can create multi page PDF's to merge to the map? Is there a reportlab method that can simply copy the receptor table to a PDF (even if it is multi page)


texts: Is ArcGIS Online unable to change the comma separator from comma to point?

texts:
I am building an ArcGIS-Online Dashboard for a German partner. I have installed indicators which has a hosted feature as it's data source. When showing the values, the decimal separator is a point which is wrong. In Germany, the comma separator is used as a decimal point, and the point is used as a thousand separator.
I can see no way of changing it. Is this problem solvable?


texts: What is causing the error "more than one linestring in KML"

texts:
I exported a KML from Mapinfo as a single combined line, however when i try to uploaded into a different mapper I get an error indicating that there are multiple linestring tags in my KML.  I can see it as a single line in google earth, and when i export it from Mapinfo its exported as a single object.

It was originally a route created from multiple segments which i combined into a single segment using the combine tool in mapinfo.

Why does it say that i have multiple line string tags?


texts:
I have one shapefile with XY coordinates in meters. Now, I want to develop one script in Python which converts meters to decimal degrees. Please help me to develop this.

Thanks


texts: Can I use my SearchCursor in Select by Location

texts:
I am trying to select the largest polygon that contains each feature in a feature class but none of the features in another feature class. I have written the code that seems logical to me but (of course) it doesn't work.

I am using a search cursor to iterate through my first feature class and for each feature I want to run my point in polygon query. it seems that my cursor is not a valid input. Is there a special way to express it, or is it not possible to use the search cursor as an input to this. If not is there an other was of searching on each point one by one?

Polygons = r"\\...\\tests.gdb\\out_contour_g"
peakPoints = r"\\...\\tests.gdb\\inv_peak"
troughPoints = r"\\...\\tests.gdb\\peaktest"
OutputFC = r"\\...\\tests.gdb\\OUTPUT1"

arcpy.RepairGeometry_management(Polygons)
arcpy.AddField_management(Polygons,"area_calc","DOUBLE")
arcpy.CalculateField_management(Polygons, "area_calc","!Shape!.area","PYTHON_9.3")

points = arcpy.da.SearchCursor(peakPoints, "SHAPE@")
for point in points:
    arcpy.SelectLayerByLocation_management(Polygons, "CONTAINS", point)
    arcpy.SelectLayerByLocation_management('out_contour_g',"CONTAINS",'peaktest',"","REMOVE_FROM_SELECTION")
with arcpy.da.SearchCursor(Polygons, "SHAPE@") as rows:
    for row in rows:
        list = []
        list.append(row[0].area)
        largestArea = max(list)
        arcpy.SelectLayerByAttribute_management(Polygons, "SUBSET_SELECTION", "area_calc =" + str(largestArea)) ####the where clause is probably wrong
        arcpy.Append_management(Polygons, OutputFC, "NO_TEST")


texts: How to create line substrings with fixed lengths from a line ignoring inflection vertices in PostGIS?

texts: I want to split a line in segments of a fixed length (e.g. 10m) like indicated in the picture. Using the approach of st_segmentize how-to-buffer-a-line-with-segments-in-postgis-using-unique-segmentsize works partially fine, but counts each inflection point as end/beginning of a line, which lead to unequal segment sizes. Using the GRASS function v.split in QGIS results in the right output v.split, so I want to know if there is a similar function/ approach in PostGIS 2.5.1.


texts: QGIS 2.0.1 connecting to ESRI WMTS jpeg vs jpg

texts:
I'm currently trying to troubleshoot an issue I have connecting to ESRI WMTS services.

Esri ArcGIS Server seems to provide the format as JPEG via the rest service, but in the get capabilites document they specify JPG - which it seems QGIS and Global Mapper don't want to recognise. Would the norm amongst mapping clients be JPG, JPEG or Both?

I've included the URL to our server if anyone is keen to try it out, I've also attached an image of what I get in QGIS

http://www.topofthesouthmaps.co.nz/arcgis/rest/services/CacheAerial/MapServer/WMTS

Thanks in Advance


texts: Development of a plugin which depends on an external Python library

texts:
I'm developing a Python plugin, the target is integrating into QGIS the functionalities of a PyPI Python library called 'elevation'.

Is there an OS independent way of installing external Python libraries and making them available to the Python plugins?

I already know the different system specific procedures for installing an external Python dependency, what I'd like to do is to distribute my plugin across different environments without making people go crazy.

What I already know

Searching the web and other similar questions it seems that there isn't an OS independent solution.

Windows users

Python is packaged and distributed inside the QGIS package, so in order to install external Python libraries you have to go through the OSGeo4W Shell and use pip from there.

related links:

QGIS standalone and Python modules
Installing Python setuptools into OSGeo4W Python
How to install 3rd party python libraries for QGIS on Windows?

OS X Users

In this case QGIS uses the built-in Python which is prepackaged inside OS X and located at:

/usr/bin/python

The issue here is that the default installed pip has some limitations and requires too many privileges.

A quick workaround is to use home-brew to install Python so that your pip can run against a user-modifiable Python framework. See the details of this approach here.

After installing all the Python libraries via pip, you just need to add their directories to the PATH variable.

Alternative way for OS X users

Inside the plugin Python code you can use the pip provided by the system in order to install the required packages. Then you can make the just installed package reachable adding it to the path.

import sys
import pip
pip.main(['install','--target=/Devel/test', 'elevation'])
sys.path.append("/Devel/test")

Homemade but working solution (tested on both Linux and OS X)

I decided to go for the alternative way since it allows to keep the library update and maintenance separated from the plugin. Every time a new user installs the plugin the latest version of the elevation library is directly downloaded and installed by pip inside a subdirectory of the plugin.

self.plugin_dir = os.path.dirname(__file__)
self.elevation_dir = os.path.join(self.plugin_dir, 'elevation')

# Checking the presence of elevation library
try:
    import elevation
except:
    pip.main(['install', '--target=%s' % self.elevation_dir, 'elevation'])
    if self.elevation_dir not in sys.path:
        sys.path.append(self.elevation_dir)


texts: Verify with WKT point is in a polygon SQL

texts:
I'm using SQL Server and WKT to deal with geographic data. What I'd like to do is have something like :

SELECT pt
FROM point_table, polygon_table
where pt.point_table WITHIN ply.polygon_table

Here pt is WKT point. Ex : 'POINT(lon lat)'

and ply is WKT polygon Ex : 'POLYGON (( lon1 lat1 , lon2 lat2 , lon3 lat3 , lon4 lat4 ))'


texts: QGIS custom UTM CRS including zone number in coordinates?

texts:
Is it somehow possible in QGIS to define a custom CRS which is based on UTM but accepts Easting coordinates that include the zone number?

The point behind this is, that I've got DXF files that I would like to include in my project. All of those are defined with complete coordinates.

When I select the UTM system with the zone, they plot way too far to the East because of the zone number.

One example point in them is

3243xxxx.371,551xxxx.924

I've found this topic on the same problem
Importing UTM coords into QGIS: Handle UTM Zone

Someone mentioned EPSG:5650 there but for me that is not the correct one...

Since the DXF data is not points only I don't really know how to convert them to use a notation without the zone number...


texts: Problem looking at tables within ESRI geodatbase

texts:
I installed the most recent OSGeo4W with QGIS 1.8.0 and GDAL 1.9.2 on my Windows 7 PC.

I can load an ESRI geodatabase (.gdb file) and I can get the window
"Select vector layers to add".
When I click on any of the layers, however, I get the message "qgis.exe has stopped working"

Any idea where the problem might be?


texts: Styling/labeling specific features in QGIS?

texts:
I need to styling/labeling some specific features in a vector layer. I thought I could use the rule-based option in the style/label window with a "field=value" expression if i have only one specific feature.

Which expression I should use if I have many specific features?


texts:
I have a raster which shows the probability that a cell will be used (raster_Use). I also have a points file that shows mortality locations (points_Morts). I have overlain points_Morts onto raster_Use and extracted the probability of use for each mortality location. I also need to extract the cell ID from raster_Use for each mortality point in points_Morts. I've tried converting the raster_Use to a polygon but I've received the following error:

  ERROR 000864 Input raster: The input is not within the defined domain.
  ERROR 000863: Invalid GP data type.

I read in another thread that a solution is to change the raster values from a float to an integer. I've tried that but the resulting raster only contained 0s because the values are between 0 and 1 and thus rounding down to 0. My next thought is to create another raster that's identical to raster_Use that contains 1...n unique values and then extract those values to the points_Morts file. I can't figure out how to do so.

Does anyone have an idea of how to do so?

@Luke provided a great answer, but unfortunately it doesn't quite work. To expand on why it doesn't, identifying the column and row numbers creates a key for the points file but there still isn't a key for all cells in raster_Use. So I know which cells points_Morts falls in, however, I don't know the ID for cells in raster_Use that don't contain points from points_Morts. I tried converting raster_Use to a points object and then implemented @Luke's solution, but the Raster to Points tool creates 4 points for every raster cell instead of 1. Ideally, I would like to start with a cell (e.g., the lower left corner cell) and assign it a value of 1 then continue onto the next cell assigning a value in ascending order.

I misunderstood Luke's question on "why?" I'm doing this. I am validating the predictive ability of a model. Raster_Use is my prediction and points_morts is my validation dataset. I'm validating my model by grouping the cells from raster_Use into 20 equal-area bins. Expected use is the predicted value for each grid cell. I place these cells into 20 equal-area bins and then summed by bin and standardize to 1. Observed use is the frequency of mortality locations in each grid cell, summed by bin and then standardized to 1. After that, I'll use linear regression to regress observed vs predicted use where the slope of the regression relationship between observed and expected provides an assessment of the predictive abilities of the model.


texts: Find intersection vertices of two polygons

texts:
I have defined two polygons and I can find whether they intersect or not.
I need to find the intersection vertices if they intersect or not ?

How can I do that ?

p1 = Polygon([(0,0), (1,1), (1,0)])
p2 = Polygon([(0,1), (1,0), (1,1)])
print(p1.intersects(p2))

This gives whether the polygons intersect but how to find the intersection vertex. I tried the following but it gave an error:

print(p1.intersects(p2).coords)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'bool' object has no attribute 'coords'


texts: Cloud masking the Landsat-8 Surface Reflectance product without using GEE

texts:
I have a collection of Landsat8 Collection2 Level2 images (surface reflectance) and I would like to cloud mask those images somehow, so I can further make a time series of NDVI and other remote sensor indices like VCI, TCI, VHI using ArcGIS pro and model builder. I tried using the Cloud masking plugin in QGIS but it doesn't load any images after I insert the MTL.txt it shows this error:

From this I understand that I am supposed to use this plugin with Level 1 product or am I driving insane? I have the same images from Level 1 and the plugin works with them. So first, is there any other way to remove clouds from the surface reflectance product using my collection? I tried some codes I found here and using Google Earth Engine but I am not very familiar with that and I am in a rush (deadline soon) so I don't have very much time to spent in GEE. And I would possibly love to find another way.  And second is it maybe possible to Cloud mask the level 1 images and then apply that mask to the surface reflectance product?


texts: Visualizing variability in the size and shape of polygons

texts:
I am working with a data set of area restrictions imposed by criminal courts. These restrictions ban a person from entering a demarcated area. These restrictions vary in their size, shape, and location. I would like to describe and, ideally, visualize variability in the size and shape of these restrictions.

My hypothesis is that area restrictions for so called inchoate offenses (public order and drug related offenses, so called "victimless crime") will have much greater variability in their size and shape compared to other offenses.

Would it be possible using PostgreSQL to modify all restrictions associated with a particular offense so that they would all share a common centroid? This would allow a visualization of the effect I am trying to describe.


texts:
Assuming you have 20+ Qgis projects each with dozens of layers:

Background: Web Map Tiles,
Features: Millions of Polygons & Points
Features are labeled, labels are never overlapping with each other
Styled to look perfect in any zoom Layer
Gives perfect results if printed out via the composer as atlas

Problem: If you want to convert the whole project with right this styling into MBTILES or Slippy Maps, what is the best solution ?

I spend for 4days testing and failing with:

QTiles (Problems with Labels, which occur in every Tile, slow)
QMetatiles (Crashes very often while saving directly to mbtiles, slow)
Converting an Atlas as GeoTIFF via GDAL2MBT to MBTILES (unstable, labouros if you have many projects)
Setting up a QGIS WMS Server and reading the tiles with the Mobile Atlas Creator (labourous if you have many projects and unstable)

Is there a universal, robust solution to export the QGIS Canvas to Tiles, fast (like a "Save As ..." Button)?

After reading every google result twice, testing many github tools, I'm blind or searched the wrong keywords. I'm open for obvious and strange suggestions, as long as it doesn't involve a reprogramming for every export.


texts: Decision tree classifying of satellite image in R?

texts:
I want to classify satellite image in R using RWeka classifier, J48. I have a CSV file with the classes required, and raster data loaded in R. I am able to make the tree, however, I am not able to use the same tree to classify my satellite image.

Here's how it goes.

inputfile=read.csv("E://R/Trial/naini.CSV")
m1 <- J48(Class~B1+B2+B3+B4+B5+B7, data = inputfile)
if(require("party", quietly = TRUE)) plot(m1)
raster<- brick("merge.tif")
plotRGB(raster, 4, 3, 2, stretch="hist")

And then I am attempting to use the tree generated to classify my satellite image, but I keep getting errors.

predict (raster, J48, filename=out_classify, progress='text',format='GTiff', datatype='FLT4S', type='response', overwrite=TRUE)
Error in trim(filename) : error in evaluating the argument 'x' in selecting a method for function 'trim': Error: object 'out_classify' not found

And if I change the arguments, this is what I get. I am not sure how to proceed from this point forward.

predict (J48, raster, filename=out_classify, progress='text',format='GTiff', datatype='FLT4S', type='response', overwrite=TRUE)
Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "c('R_Weka_classifier_interface', 'R_Weka_interface')"

I am new to R.


texts:
How do I create such map from point data?

I am using ArcGIS 9.3 and 10


texts: Is there a method in PHP to covert a shapefile to GeoJSON? I need to upload shapefile and load it into a PostgreSQL table. I tried PHP exec(), but it didn't work.


texts: Logically assigning values to polygons starting at top left using ArcGIS Desktop?

texts:
I have a layer in which the order of the OBJECTID's is rather random.

How do I assign a new field in the data table that assign's number to each polygon based on it's position on the map?

I'd like to have the top left polygon have number 1 and then increment from left to right, top to bottom.


texts: Can spatial data in Oracle be imported into SQL server?

texts:
Oracle supports spatial data and so does SQL Server.

What are the things to consider if we have to migrate spatial data from Oracle to SQL Server?
Are there any features in Oracle that are not there in SQL Server for supporting spatial data?

We're talking about data maintained by ESRI ArcGIS if that affect answer of the question.


texts: Setting transparency for all layers in a data frame in ArcMap using ArcPy?

texts:
Using ArcMap 10.6, I've tried the following code from a previous post and get an error or it executes but does nothing. Here is the code and error I receive:

Code

mxd = arcpy.mapping.MapDocument ("CURRENT")
df = arcpy.mapping.ListDataFrames (mxd, "Layers") [0]
for lyr in arcpy.mapping.ListLayers (mxd,"", df):
    lyr.transparency = 0
arcpy.RefreshActiveView()
mxd.save
del mxd

Error

  RuntimeError: LayerObject: Set attribute transparency does not exist


texts: Cloud mask error: Pattern 'cfmask' did not match any bands

texts:
I want to "replace" the values of the cloud pixels for: it could be a mean, to calculate pre and pos-date NDVI for one country.
I've used the code of a previous answer:
Mask clouds in LandSat 8 surface refletance image

But this is giving error:

ImageCollection (Error)
Error in map(ID=LC08_202033_20170804):
Image.select: Pattern 'cfmask' did not match any bands.

What am I doing wrong?

The code:

//Choose country using GEE Feature Collection

var region = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw').filterMetadata('Country', 'equals', 'Portugal');

//Add region outline to layer ‐ for selected countries

Map.addLayer(region,{}, 'Portugal');

// image collection pre 11/10/2017

var lt8_pre = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filterBounds(region)
  .filterDate(ee.Date("2017-08-01"),ee.Date("2017-10-11"));

var lt8_ndvi_pre = lt8_pre
  .map(function(img){
  return img.addBands(img.normalizedDifference(['B4', 'B5'])).updateMask(img.select(['cfmask']).neq(4));
  });

print(lt8_ndvi_pre);

Map.centerObject(region, 10);
var ndvi_viz = {min:-0.8, max:0.3, palette:'000000,00FF00'};
Map.addLayer(lt8_ndvi_pre.select('nd').mean(), ndvi_viz, "LT8 2017 NDVI mean pre 11/10");


texts: How to force extra dimensions on MultiLineString geometry type?

texts:
I attempted to use ST_Force4D() I've used for MultiPolygons (and other geometry), but it doesn't seem to be the correct function.
INSERT INTO table_one (geom) SELECT ST_Force4D(geom) FROM table_two;

HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
********** Error **********
ERROR: function st_force4d(geometry) does not exist

How can I force geometry of type MultiLineString to type MultiLineStringZM?


texts: Do local Bing Maps servers or Google Earth servers exist?

texts:
Is it possible to have all the map data come from local server instead of Google Earth or Bing Map but still use their tools as front end.

For example Bing Map control (from their SDK) would be talking to some local server in the network. I want this kind of setup in a place where there is no internet and the client wants to use Bing Map toolset. Same question applies to Google Earth.

OR is it that its only possible with something like ArcGIS like server software but in that case Google Earth or Bing Map control won't be its front end.


texts: Why QGIS 2 and QGIS 3 return different value for $area function?

texts:
I measure the area of a polygon with the QGIS Calculator: The cartesian area($ geometry) calculation function gives the same result for QGIS versions 2.16 and 3.4.5, but it is different with the $area function as shown in the image:

I don't understand why because the EPSG code and the ellipsoid are identical.

Could someone explain to me?


texts: Populating oracle sdo_geometry object using GeoJSON?

texts:
I have GeoJson data. I would like to convert this  GeoJson data into sdo_geometry object type and insert Oracle spatial database. Which columns represent geojson data  in oracle sdo_geometry object.

For example,
My GeoJson Data is

  {"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Point","coordinates":[410694.5082851944,4466102.355016577]},"properties":null}]}

sdo_geometry in Oracle is consisted of columns that are named SDO_GTYPE, GEOLOC.SDO_SRID, GEOLOC.SDO_POINT, SDO_ELEM_INFO,SDO_ORDINATES.

How Could I populate columns which belong to sdo_geometry object in Oracle using GeoJson object.

Geojson objects might be different, so it couldn't always point type and because of this, geojson object could have different values.


texts:
I have a Alberta province (Canada) section polygon shapefile in GCS_NAD83 with 255,000 polygons. I like to split each section into 4 equal quadrants (NW, NE, SW, SE)

The Alberta province has several UTM zones, so I guess I can convert it to Lambert Conformal Conic or Albers Equal Area CS if the section scan be split only in projected units.
There are people selling this information, but I was wondering if this is something I can create using the data I have?
Please let me know. Thank you.


texts:
I'm trying to create points on the edges of polygons.

So far I have converted the polygons to lines, then exploded the lines, then ran the following code which plots a point in the middle of the lines.

layer = iface.activeLayer()
temp = QgsVectorLayer("Point?crs=epsg:27700", "result", "memory")

temp.startEditing()
attrs = layer.dataProvider().fields().toList()
temp_prov = temp.dataProvider()
temp_prov.addAttributes(attrs)
temp.updateFields()

for elem in layer.getFeatures():
    feat = QgsFeature()
    geom = elem.geometry().interpolate(elem.geometry().length()/2)
    feat.setGeometry(geom)
    feat.setAttributes(elem.attributes())
    temp.addFeatures([feat])
    temp.updateExtents()

temp.commitChanges()
QgsProject.instance().addMapLayer(temp)

I have some polygons which are adjoined, sometimes several squares in a row.

I've noticed that each polygon still has a separate line, even though they're adjoining polygons if that makes sense?

I want to remove the points where two polygons adjoin, and so therefore I just want to remove all instances where there's duplicate points/lines. So not just remove duplicates so there's one left, but remove both instances of duplicates.

Hope that makes sense and someone can help!


texts: Check whether QGIS3 was started by qgis.bat or qgis.exe

texts:
Is it possible to check whether QGIS3 was started by executing qgis.bat or qgis.exe directly?

I want to do some checks in the qgis.bat before qgis starts and make sure that qgis is always started by calling the qgis.bat.

Edit1:

One idea is to set an environment variable in the qgis bat and check in the startup.py if this environmental variable is found:

bat-file:

set batchfilestart="YES"

python part

import os
try:
    env_value = os.environ.get("batchfilestart")
    if env_value !="YES":
        print ("QGIS not started by calling the bat file")
        print (env_value)
    else:
        print("QGIS started by calling the bat file")
        print (env_value)
except KeyError:
    print("Not exist environment value for %s" % "key_maybe_not_exist")

Is there a better way to make sure I can always run some checks before QGIS starts?


texts: ogrListLayers() in rgdal returns character(0) for a kml file

texts:
I try to read a kml file in R. In the past, there was a post (link). I checked stack overflow and found similar posts. By following the examples, I was trying to read a kml file. I could read the kml file using getKMLcoordinates() in the maptool package and extract geo data (i.e., longitude and latitude). I also used the rgdal package in order to see layer names. But, ogrListLayers(dsn = "name.kml") surprisingly returned character(0). My questions are:

Is is possible for kml files not to have layer names?
If there is no layer name, is there still any hope to find the names? Or do you have to simply give up?

P.S. I originally posted this message in stack overflow, but I have not received any suggestions. So I decided to post this message here.


texts: Specifics of ESRI's use of FlexNet Publisher TrustedStorage at ArcGIS 10

texts:
So at ArcGIS 10, ESRI adopted Flexera's FlexNet Publisher Trusted Storage module for secure exchange of license file "fulfillment records" and transfer of product activations.

Question what network traffic occurs and what corresponding firewall ports should be open to allow the traffic to pass? Also, what conflicts will arrise with other applications running the same FlexNet Publisher licensing tools--Autodesk, SolidWorks, Adobe, etc. and how best to make things coexist?

The FLexNet Publisher License Administration Guide doesn't address those specifics. And when looking at the ArcGIS 10 installations you'll now see a FLEXNetLicensingService.exe running and calls to libFNP.dll on license server and client. That in addition to the "old" FlexNet Publisher lmgrd.exe and ESRI's FlexNet ARCGIS.EXE Vendor daemon on the license server.

How does it all behave? And how do we fix it when collisions or version conflicts occur?

A very scary thing is that a number of Apps using FNP are placing configurations and executable files into the common directory:

C:\Program Files\Common Files\Macrovision Shared\FLEXnet Publisher

Where vendors can do what they want and update their FlexNet components at their convenience, without regard for release level of other applications present on the system. Looking at the FLEXERA site I don't have a sense that they make much of an effort to support the end users. License exchanges over VPN has been especially troublesome for end users.


texts: Adding data to Web AppBuilder in ArcGIS Online?

texts:
I'm trying to use Web AppBuilder in ArcGIS Online to create web application. I can start the Web AppBuilder and everything looks fine. But when I try to add data to the application from my ArcGIS Server or from my existing  webmap on ArcGIS Online, I can't see any data. It looks like the Web AppBuilder is trying to add data but nothing happens?

Any ideas ?


texts: Calculate total coverage of circles on map with PostGIS

texts:
I am trying to calculate the total area covered by c.a. 400 Circles of diameter 15Km. The centers of these circles are stored in a column of type geometry:4326 in my database.

Some circles overlap, and I should not add to the total area the overlapping amount twice.

Until now I tried the follwing approach:
1) Use buffer to get circles: select st_buffer(get_location, 15000) from d_station;
2) Use st_union to join the circles into a single multipoint geometry
3) use st_area to calculate area.

Unfortunately this approach doesn't work, I think because I am not using st_union in the correct way?

Then I found an interesting tutorial here, which explains something very similar to what I need. Unfortunately I am not able to get a useful result using this approach either.

How can I obtain a list of circles and merged circles polygons, and then calculate the sum of their areas?

Update: I am trying Pierre suggestion. Here the steps I performed:

CREATE temporary TABLE pts AS
select ds.id, geo_location as geom from pa.d_station ds

CREATE temporary TABLE circles AS
SELECT id, ST_Buffer(geom, 15000) AS geom FROM pts;

create temporary table circles2 as
SELECT a.id, pa.ST_DifferenceAgg(a.geogr, b.geogr) geogr
FROM circles a,
     circles b
WHERE ST_Equals(a.geogr, b.geogr) OR
   ((ST_Contains(a.geogr, b.geogr) OR
    ST_Contains(b.geogr, a.geogr) OR
    ST_Overlaps(a.geogr, b.geogr)) AND
   (ST_Area(a.geogr) < ST_Area(b.geogr) OR
    (ST_Area(a.geogr) = ST_Area(b.geogr) AND
     ST_AsText(a.geogr) < ST_AsText(b.geogr))))
GROUP BY a.id;

select sum(st_area(c2.geogr)) from circles2 c2;

But I get a result which is way too low from the expected. Is my approach wrong?


texts:
I have a QGIS project that makes heavy use of WFS for source layers. I tend to prefer this because I like to ensure that my map will stay current with changes to authoritive data sources (e.g. roads, administrative boundaries). My current region of interest is very small (a single town and small hinterland), but the WFS coming from a national government source is very large (many MBs of features that I'll never look at). QGIS seems good at limiting my requests, but on loading my project I find that it has to spend several minutes making WFS requests and re-downloading data that I'm positive is not going to change in the period that I'm working on my project.

I find this time waiting on these large WFS queries to complete is frustrating. Occasionally waiting on the response actually causes QGIS to crash, and I have to re-open QGIS and then it will start the downloads again from scratch.

Is there a way to set explicit client-side caching for WFS layers in the client (QGIS) to prevent a re-download? If I could tell QGIS to only perform a new query after a month (or whenever I deliberately clear or invalidate the cache), that would be ideal. I want this to work across sessions.

Existing similar questions and why they're not relevant:

WFS Feature Caching in QGIS | Actually about a cache feature limit problem in the same QGIS session rather than caching WFS in QGIS across sessions.
WFS Caching / Back button | Caching in the browser.

I'm using QGIS 2.99 (master), which I gather is more or less QGIS 2.18.

Something I just observed: I found the reloading process this morning to be much faster than last night when I wrote a draft of this question, so it could be that I'm mistaken and that there already exists some kind of client-side cache or perhaps a mechanism for producing a diff and making a targeted query. However there was still one layer (landuse: with many thousands of complex features) that appeared to completely re-download. Is this then perhaps a matter of missing server-side cache headers for this particular service? This layer is from a different source than the rest of my layers. But is there still a way to declare a cache expiry in QGIS as the client, even if the server doesn't adequately set cache headers?

Indeed this question (Cached Featured in WFS - can I change default) suggests that QGIS does indeed default to caching WFS responses. (But the question asks whether there is a setting to override the use of a cache, which is related... and unanswered).

These questions (Is the QGIS cache directory automatically created? and How do I use a filter on an uncached WFS layer in QGIS?) tell me that QGIS definitely does caching. I'm only surprised because I keep seeing (some) of my own layers re-download after I restart.


texts:
I have an ECEF position as:

x: 1776694.9644800487,
y: 5759463.252896432,
z: 2079310.7706825584

Using some online ECEF to LLH converters, this value is transformed to:

19.152548, 72.855909, 8.00

I am curious how is 2079310.7706825584 transformed to 8.00.


texts: Problems with specifiying TM_fill when creating a map in R with TM_maps

texts:
For my research I am visualizing plastic contents in segments of water.
After creating my segment (online, polygons in KML-file) I jointed it with my table with fieldworkdata with plastics per segment.
Creating a map with these segments works fine until I want to combine the amount of plastic and the segments. My code is the following:
rm(list=ls())
setwd("C:/Users/......................../Thesis") # Set working directory (change to your own location)

#packages
library("sf")
library("raster")
library("sp")
library("spData")
library("spDataLarge")
library(tmap)    # for static and interactive maps
library(leaflet) # for interactive maps
library(ggplot2) # tidyverse data visualization package
install.packages("rgdal") #importing KML files
library("rgdal")

(rgdal_show_exportToProj4_warnings="none")

library(dplyr)
library(stringr) # for working with strings (pattern matching)
library(tidyr)

#Load data for segments
Meet3 <- read.csv2("Weeklydata/DataMeet3.csv", header = TRUE)

segmenttest <- st_read("R maps/polygons/PolygonsLeiden.kml")
Meet3$Segments <- c("1","2","3","4", "5", "6","7","8","9","10", "11", "12","13","14","15","16", "17", "18","19","20","21","22", "23", "24","25","26","27","28", "29","0")
PlasticData =Meet3[c(-30),]

PolOnly <- segmenttest[ ,1,3]
PolOnly$Segments <- c("1","2","3","4", "5", "6","7","8","9","10", "11", "12","13","14","15","16", "17", "18","19","20","21","22", "23", "24","25","26","27","28", "29")

#Combine data
PlasticMap = left_join(PlasticData, PolOnly)
current.mode <- tmap_mode("view") #switch to interactive map

plot(st_geometry(PlasticMap$geometry), col = PlasticMap$total.per.segment) (this works fine!!, but i want to create TM_maps)

Map1 = tm_shape(PlasticMap$geometry) +   tm_polygons("total.per.segment")

The problem is that when running Map1, the following error occurs:

Error: Fill argument neither colors nor valid variable name

Link to 2 files with example of the datafiles:
https://we.tl/t-jLvLZHlWUM


texts: Geocoding with MMQGIS creates empty attribute table from CSV file

texts:
I need to geocode an Excel file of addresses in Greece ("number", "name", "address", "city", "country" are my columns) using MMQGIS with OpenStreetMap.
I took all the steps to turn it to a .csv-file (converted .xls to .txt, then to UTF-8 .txt, then to .csv). I managed to create valid locations to save my output files, ran the Geocoding, but, my shapefile -as well as the attribute table- is empty of content.

Any ideas and suggestions would be very welcome!
PS: My actual Excel file is quite longer, this is a file for practicing.


texts:
The tool I'm developing is written in Java with Spring Boot, using PostgreSQL as database (PostGIS extension). I also make use of GeoTools.

I want to calculate the geometrical difference between two LineString geometries, but I'm quite lost in this problem.

Perhaps it'd be better to draw an example.

The red LineString would be Trip 1:

The green LineString would be Trip 2:

What I want to get is the orange LineStrings (as MultiLineString or so):

The problem is that to get the difference between red and green lines I'd like to apply some kind of tolerance values, e.g. 5 meters, so that I can get rid of inaccuracies like this one (and consider them as lines that crosses the same path):

In PostGIS, there is a function called ST_SymDIfference, but there isn't tolerance parameter.

Is it possible to accomplish what I want?

Edit

@eth4io answer is the correct one. Since its based on PostGIS code, I'll write the Java code which replicates it:

public Geometry calculateDiffs(final LineString pPlan, final LineString pRoute) {
    Double distance = 0.0002d;
    Geometry planBuffered = pPlan.buffer(distance);
    Geometry routeBuffered = pRoute.buffer(distance);
    Geometry intersection = routeBuffered.intersection(planBuffered);
    Geometry difference = pRoute.difference(intersection);
    return difference;
}


texts: Get custom records from Address locator service in ArcGIS Server

texts:
I've already created an Address Locator successfully and published it as a Geo coding service, but have some additional requirement.

Let's say I've a field in my data set named category which have different values like

Houses
Parks
ATM

What I want is to pass this category value as parameter and get the service results from only those records which match the category parameter value.

I need to know how to make it happen?


texts: Find minimum elevation along a polyline and return xy or point?

texts:
Using ArcGIS 10.1 and Python IDLE 2.7.2...

I'm trying to find the minimum value of a polyline with a raster layer and turn that LOCATION into a point. I put everything together in model builder, but the raster calculator doesn't like 32-bit floating points, which is what my DEM is. I tried to get around the raster calculator with the Con tool in python, but I believe an attribute table is required for this (which I don't have due to floating point). I could convert the DEM to integer and build one, but this is already a lengthy process that involves many scripts and tools and wanted to see if anyone else out there had any thoughts? I'm open to any anything, really, conversion of polyline to polygon, etc...attaching a screenshot of modelbuilder process that works, but I need a script as I am doing this 1000+ times.

I can turn the polyline into a 3D feature and find the minimum value, but I don't know how to find WHERE that minimum value is. If I could locate it with XY coordinates, I can then generate the point.

Thoughts? If I'm wrong about the Con tool and my syntax is wrong, that would be awesome, so please feel free to suggest something. It finds where the min is along the polyline that has been converted to raster and where that min is within the DEM, thus returning a min pixel:

tempRas1 = Con(Raster(LG_DEM_0_tif == LG_0_Zonalmin_img), LG_DEM_0_tif)


texts:
Is there a way to use ArcPy without having ArcGIS Pro?
I am currently using ArcMap 10.8 and am trying to set up a conda environment that is compatible. However, I'm having trouble finding the esri packages I need to install and all tutorials I find are made specifically for ArcGIS Pro.
I found a couple solutions here: Making separate Python installation that can call ArcPy?,
but they're pretty dated. I am considering trying the provided script:
try:
    import archook #The module which locates arcgis
    archook.get_arcpy()
    import arcpy
except ImportError:
    # do whatever you do if arcpy isnt there.

but if there is an esri-approved way to import arcpy without Pro, I'd prefer that. I just can't find it.


texts: Labeling only points that ends with the number 5 in QGIS

texts:
I have a shape file with points at evenly spaced intervals, the table has the following fields:
Name  is up to 32 characters of text.  it is currently populated with an eight digit number.

Point_No is an integer field which has the same eight digit number that is in the Name field.

What I would like to do is to only label the points that have the Name field or the Point_No field  which end with the number "5"

I have tried case statements, substr statements, I have tried to use the modulus with the Point_No field, but I keep getting parsing errors.


texts: QGIS Python Error : No module named 'encodings'

texts:
I have written a python code which uses the 'Raster Surface Volume' algorithm of QGIS Desktop 3.8.3. When I execute the python code, I am getting an error :

The system cannot find the path specified.
Fatal Python error: initfsencoding: unable to load the file system codec
ModuleNotFoundError: No module named 'encodings'

Current thread 0x00002218 (most recent call first):
Press any key to continue . . .

I have imported the libraries as given below.

from osgeo import ogr, gdal
from gdalconst import *
from qgis.core import *
import qgis.utils, sys
from qgis.gui import *
from PyQt5.QtGui import *
from PyQt5.QtCore import Qt
from PyQt5.QtCore import QFileInfo
app = QApplication([])
QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
QgsApplication.initQgis()
sys.path.append('C:/OSGeo4W64/apps/qgis/python/plugins')
import processing
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *
import os, glob

In the 'Environment Variables', in the 'path' of 'user variables' section, I can find :

C:\Users\sreeraj\AppData\Local\Programs\Python\Python37\Scripts\
C:\Users\sreeraj\AppData\Local\Programs\Python\Python37\

I can also find the same in the 'path' of 'system variables' section in the 'Environment Variables'.

In the folder 'C:\OSGeo4W64\apps\', I can find 2 folders namely 'Python37' and 'Python27'.

I have installed 'python-3.7.4-amd64.exe' in my computer before installing the 'osgeo4w64'. Which means that I am using Python 3.7.4 in my computer and when I installed 'OSGeo4W64', I can find find a folder 'Python37' (along with 'Python27') which is integrated to QGIS Desktop 3.8.3.

Can anyone help me with the exact solution to this error 'No module named encodings' ?

As @user2856 suggested, I opened OSGeo4W command prompt, then executed py3_env, qt5_env and ran the script using python3 path\to\your\script.py . Now, I got a QGIS import error :

from qgis.core import * ModuleNotFoundError:
No module named 'qgis' .

python-qgis.bat file in C:\OSGeo4W64\bin\ is given below :

@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 QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%PYTHONPATH%

"%PYTHONHOME%\python" %*

I have tried the accepted answer of @Germán Carrillo of Problem with import qgis.core when writing a stand-alone PyQGIS script by running var_entorno_osgeo4w.bat in OSGeo4W64 shell as given below :

run o-help for a list of available commands
C:\OSGeo4W64>var_entorno_osgeo4w.bat
C:\OSGeo4W64>set PYTHONPATH=C:\OSGeo4W64\apps\qgis\python
C:\OSGeo4W64>set PATH=C:\OSGeo4W64\apps\qgis\bin;C:\OSGEO4~1\apps\Python37\Scripts;C:\OSGEO4~1\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\WINDOWS\system32\WBem
C:\OSGeo4W64>

Still, I am unable to solve this issue. How can I solve this QGIS import error : No module named 'qgis' ?


texts: How to set "Spatial ref" for table in GeoPackage

texts:
I have a GeoPackage table with no crs definition. I would like to fix this now.
There is no command like <ogrlayer> .SetSpatialRef (...)

Download d.gpkg

# -*- coding: utf-8 -*-
import osgeo.ogr as ogr

driver = ogr.GetDriverByName("GPKG")
data_source = driver.Open("c:/..../d.gpkg", update=True)

ly_25833=data_source.GetLayerByName("tab_25833")
ly_without=data_source.GetLayerByName("tab_without")

sr_25833=ly_25833.GetSpatialRef()     # EPSG:25833
sr_without=ly_without.GetSpatialRef() # None

print (sr_25833)
print (sr_without)

# !!! what i wont !!!!
# sr_without=ly_without.SetSpatialRef(sr_25833)

data_source.Destroy()

Download d.gpkg


texts: find latlngs that are some x distance away from a point to south to north to east to west

texts:
I have a latlng i want to find four points with equal distance from this latlangs in all directions east, west, south an north is there any formula .

I am trying to build polygons around higways so draw polygon around poly line the distance will be less than 1 km.


texts: Turn off field programmatically in ArcMap via VB.Net

texts:
I'm creating an addin button in Arcmap via VB.net that when clicked, runs the identify tool. But I want only selected fields to appear in the identify dialog. My idea is to hide the unnecessary fields and then show them again once the identify tool is deactivated.

Is it possible to programmatically turn off/on a field?


texts: Lat and long extents of Canadian provinces

texts:
I'm seeking the latitude and longitude extents (i.e., extremes of lat and long) for each Canadian province. I already have Alberta and BC.  For example, Wikipedia lists Ontario's extent as 49° 00’ N, 114° 04’ W to 60° 00’ N, 139° 03’ W.

I need the same for the other provinces. I googled a variety of terms.


texts:
I am trying to get a kml from a URL to display in OpenLayers. However, it is a ground overlay so it cannot be displayed the standard way. I found this thread:
Does OpenLayers supports KML Ground Overlays?
However, I am very lost as to how to parse this manually. I know I can make an AJA request and parse it, but don't think I can do that bc of cross-domain issues.

So could someone guide me in the right diirection to get a KML parsed directly from a url?

Thanks!


texts: Disable QGIS Tips in Python on startup file

texts:
I want to create a portable QGIS. It can be done by copy the folder of QGIS2 and put in other computer. And if I want to run python in startup, I have to put the startup file in C:\Users\USER\.qgis2\python.

So, I can't do it the way in : Disable QGis Tips Panel when application start

I want to do this in python. I have to input the script in startup file. So, how could I script this?


texts: Error when trying to create a DOY band on an ImageCollection in Google Earth Engine

texts:
I posted this in the developers forum as it sounds to me like an internal error rather than an error in the code but I've not received a response so I'm trying my luck here.

I'm trying to create a DOY band on an image collection. I found a similar question on stackexchange which includes a function for doing this (Add a date (day of year) band to each image in a collection using Google Earth Engine). This code works fine.

When I try to run the function on my ImageCollection I get an error message:

ImageCollection (Error)
Error in map(ID=20170601T102019_20170601T102022_T33UUU):
Image.date: Image does not have a valid system:time_start property!

This sounds to me like there is an error in the metadata in one of the images in my collection. Can anyone help solve this?

NB. 'aoi' is a multipolygon shapefile.

Here is a link to my code in gee:

https://code.earthengine.google.com/f4fc313d1fa22f5ff6e0a9eb802f8970


texts: Using Mapbox Studio (aka Tilemill 2) to create a self-hosted map

texts:
What would be the neccessary prerequisites to run a map/tile-server (serving raster tiles) using an exported vector source and tm2z-styleset?

Has anyone set up a working environment for that?

Is Tessera a possible package for such a case?


texts: Converting part of raster into .xyz file

texts:
I have a very large raster whose GDAL Info is as below:

Size is 172800, 67200
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,83.999167206000010)
Pixel Size = (0.002083333000001,-0.002083333000001)
Metadata:
AREA_OR_POINT=Area
ATTRIBUTE_LABEL=BDTICM_M
ATTRIBUTE_MEASUREMENT_RESOLUTION=1
ATTRIBUTE_TITLE=Absolute depth to bedrock (in cm)
ATTRIBUTE_UNITS_OF_MEASURE=cm
CITATION_ADDRESS=shanggv@bnu.edu.cn / tom.hengl@isric.org
CITATION_ORIGINATOR=College of Global Change and Earth System Science, Beijing Normal University / ISRIC - World Soil Information
CITATION_URL=http://dx.doi.org/10.1002/2016MS000686
CONFIDENCE_INTERVAL=M
DATA_FORMAT=Byte
DEPTH=None
DOWNLOAD_FTP_URL=ftp://ftp.soilgrids.org/data/recent/BDTICM_M_250m_ll.tif
HORIZON_LOWER_DEPTH=5000.00 m
HORIZON_UPPER_DEPTH=0.00 m
KEYWORD1=bedrock
KEYWORD2=depth
NO_DATA=255
PROJECT_URL=http://soilgrids.org
PUBLICATION_DATE=2017-03-10
RANGE_DOMAIN_MAXIMUM=5000
RANGE_DOMAIN_MINIMUM=0
SERIES_NAME=SoilGrids250m
SLD_NAME=soilgrids250m:BDTICM
SLD_URL=ftp://ftp.soilgrids.org/legends/BDTICM.sld
TECHNICAL_SPECIFICATIONS_URL=https://github.com/ISRICWorldSoil/SoilGrids250m
Image Structure Metadata:
COMPRESSION=DEFLATE
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-180.0000000, 83.9991672) (180d 0' 0.00"W, 83d59'57.00"N)
Lower Left (-180.0000000, -56.0008104) (180d 0' 0.00"W, 56d 0' 2.92"S)
Upper Right ( 179.9999424, 83.9991672) (179d59'59.79"E, 83d59'57.00"N)
Lower Right ( 179.9999424, -56.0008104) (179d59'59.79"E, 56d 0' 2.92"S)
Center ( -0.0000288, 13.9991784) ( 0d 0' 0.10"W, 13d59'57.04"N)
Band 1 Block=172800x1 Type=Int32, ColorInterp=Gray
NoData Value=-32768
Overviews: 86400x33600, 43200x16800, 21600x8400, 10800x4200, 5400x2100, 2700x1050, 1350x525`
I want to extract from it a part of data as xyz file. But I cannot seem to know how to do it. I have tried `gdal_translate -of XYZ -a_ullr -89.2972 60.0378 -117.719 44.901 "D:\\Mesh\\Shangguan et al. (2016) Bedrock data\\BDTICM_M_250m_ll.tif" "D:/Mesh/MRB/ascii shangguan for NCRB 0125 .xyz"

but it writes an xyz file more than 100GB and the size is increasing even more than 100GB as I write the question. Can I just have the data for "-89.2972 60.0378 to -117.719 44.901"?


texts:
I've looked all over the net and I'm wondering if there are any plans by the QGIS developers to integrate some of the newer versions of SAGA available in the repository and as binaries on the project's homepage into QGIS 3?

In a few instances, I've run into issues where an algorithm available in a more recent version of SAGA merits me installing it simply so I can have access to its tools. Having two installations of SAGA is a counterintuitive though.

Can anyone point me to a discussion thread that has possibly come up?


texts: Automatically build branches of a river network?

texts:
A small non-profit specializing in fish habitat restoration, has asked me to link a spreadsheet of field notes to their GIS.

What we have is:

A spreadsheet of field notes for small fish streams.
A shapefile of unnamed rivers. Each branch of a river is often
comprised of several individual lines.

What we need is:

Merge together lines that form a single branch of river.
Give a name to every branch of river corresponding to the spreadsheet
notes.

Ideally, I will pan across the shapefile, clicking each branch and entering data for the "Name" attribute. I do not want to have to mess about, zooming in and merging multiple lines into a single branch, before I name it.

Using QGIS or ArcGIS extension or script.


texts: Extracting the standard error map from autoKrige in R

texts:
autoKrige in R works fine and plotting the output gives a good overview of the interpolation with standard error. Following these steps:

library(sp)
library(automap)

loadMeuse()
demo(meuse)
# Ordinary kriging
kr = autoKrige(zinc~1, meuse, meuse.grid)
plot(kr)

#Extracting parts from the autoKrige object
prediction_spdf = kr$krige_output
sample_variogram = kr$exp_var
variogram_model = kr$var_model
plot(prediction_spdf)

In the last line I only plot the extracted prediction. I also want a separate standard error plot, but I can't extract that from the autoKrige output. Anyone who knows how?


texts: Converting ASCII file to shapefile in ArcGIS Desktop?

texts:
I am trying to convert ascii  to shapefile I have huge points, latitide, longitude and elevations, I want to use geoprocessing tool also to convert this please let me know if anything is avalible.

I have text file contains longitude, latitude and elevation points I want to convert from the text file to shapefile using ArcGIS Desktop 9.3.1.

I have downloaded a tool asciitoshape but its not working properly, GUI for ascii to shape file is working, but gp.ascii2shape is not working.


texts: Match polygons of two independently created shapefiles to create dataframe

texts:
I have two independently created shapefiles

A) shapefile 1 lists the german postal codes and can be downloaded from https://blog.oraylis.de/2010/05/german-map-spatial-data-for-plz-postal-code-regions/
'plotting' will create a map with postal codes as polygons.

B) shapefile 2 (PolygonSamples.shp) contains a list of polygons which are distributed all over Germany. The created polygons might be smaller, equal or larger than the polygon of a single postal codes polygon.

How to 'match' the two shapefiles to get a dataframe that lists which polygon in shapefile 2 matches which postal codes(s) of shapefile 2?

The result would ideally look like:

 Polygon ID (shapefile 2)     Postal Code (shapefile 1)
     1                                80995
     2                                80997
     2                                80999
     3                                81247

The polygon shapefile was created in GeoMap. Not very commonly used unfortunately. I can open it in QGIS desktop, alas whatever I try in QGIS (join, etc) - saved files will always be empty.


texts:
I am part of the technical team of the Autoridad de Fiscalización y Control Social de Bosques y Tierras of the Estado Plurinacional de Bolivia , and we are conducting a multi-temporal study of the dynamics of fires in national territory, we would be very useful to access MCD64A MODIS products ; MCD64A1 , we are now mapping products and MCD45A1 MCD45A grateful,

When we try to use the tool "Modis Reprojection Tool" with MCD64A1 , the HDF data files haven't datum, and that genres error datum, so which datum we should use in HDF Products to get products geotiff with projection Geographic datum WGS84?


texts: Add locations individually in map view to existing table rows?

texts:
I have a table in which many records have long/lat info, but many others do not.  In other words the locations for those records are currently unknown.  Therefore there are many records that are visible in table view, but not in map view.  This is not an error, it is just the current state of my database.

I would like to be able to add the locations manually, case by case, as they become known. I was hoping there was a way to do that. I know you can add new records via the map view and table view. And you can edit existing records in the table view. I am wondering if you can edit an existing record in the map view by clicking on the map to indicate that location. The long and lat where clicked would populate the table.

If not, through the existing UI, what would be the best way to implement such functionality?


texts: Spatial interpolation of hourly temporal resolution temperature data using QGIS

texts:
I would like to spatially interpolate gridded (10X10km) hourly temperature data to a finer spatial grid (e.g. 5X5km). I am talking about 1 year of data (thus 8760 hourly temperatures) for each point of my grid (e.g. 100,000 data points).
I am using QGIS for the interpolation.
First, I have issues in uploading the .txt file in QGIS. The QGIS crashes as soon as I am trying to upload a .txt file with 8762 coloumns (e.g. lon, lat, T1, T2, T3, ..., T8760).
Is there a way to easily uplodad and work on heavy files in QGIS?
Secondly, I have problems with the interpolation. Using for instance IDW, what I get is a raster for each hourly tempererature time-step. On this regard, I have two questions:
1 - Is there a way to automate the procedure using QGIS? Thus, automatically obtaining hourly temperature data interpolated to my finer grid?
2 - Are there other easier procedures/software to obtain the interpolation to a finer grid?


texts: AGOL buildingexplorer. Can't show all the layers of my 3D Models at once, only 1 model at the same time

texts:
I have 2 3D models in my scene. I want to show all the layers of these models but for some reason I can only see all the layers of 1 model at the same time. I think it has to do with the building explorer tool. When I click on the tool and then on a model. It shows all the layers of this model but doesn't show all the layers of the other model.

Above: How it is normally shown

Below: When I select this building, all the layers show. But the other model has just 1 layer visible instead all of them


texts: Should I use a specific patch to upgrade ArcGis Server from 10.2 to 10.2.2

texts:
http://support.esri.com/en/downloads/patches-servicepacks/list?productid=66&productVersions=10.2.2&categoryTypes=8 Here I have fouund a list of patches and service packs. First of all whats the difference between patches and Service Packs.
Does any patch of 10.2.2 performs the upgrade, or I have to install all the patches ?


texts: Programmatically edit/update metadata in ArcGIS?

texts:
Has anyone succeeded in programmatically updating metadata in ArcGIS 10? Considering using Python/arcpy but ArcObjects (C# or Python/comtypes) is also a possibility.

I need to update both the FGDC and the ArcGIS-ISO format metadata, and whatever solution is used needs to be able to retain the existing (non-blank) elements along with the added elements, except where they are in conflict in which case the added elements overwrite the existing elements.


texts: Getting Search Widget response's lat long in ArcGIS API for Javascript 4.x?

texts: How would I retrieve the lat/long of the city/location selected using the Search Widget?  I do not need to display this location on a map.  All examples I have so far come across attach the Widget to a SceneView or Mapview.


texts: What does double digitized polygon mean?

texts: I know what is a polygon but when a polygon is said to be double digitized?  I failed to find a suitable meaning of double digitized after googling. I understand double digitized means digitized twice, so I hope it would be enough if someone explains what does digitized polygon mean. I found the term double digitized from source . ESRI Shapefile Technical Description (page 1)


texts: OpenLayers 4 add/remove features/interactions

texts:
How can I easily add/remove interaction in OpenLayers 4?

I need to add:

get layer Info as default (WMS feature on OL2)
Measure tool. (Already implemented OpenLayers4 stop interaction but I need an easy way to activate/deactivate the feature)
Plot a map marker.
Draw a circle.
Draw a polygon.

Here a demo with OL2 of what I need to do.

http://www.traffwebdemo.co.uk/rtcplus/

What I need to know if exist some kind of add feature based on which button I click and when I finish with this action (just one interaction) return automatically to the default interaction.

Example:

the default is getInfo.
click on draw a circle button
the feature is activate
I draw the circle
automatically the draw feature is removed and getInfo is turn on


texts: Creating prepared geometries with OGR from Python

texts:
Is it possible to create a prepared geometry object using the Python SWIG geos.ogr module? It seems that Shapely and GeoDjango both have this functionality, but I can't see any mention of it in the standard Python module source or documentation. Unless I'm missing something, this seems somewhat of an oversight.

http://toblerity.github.io/shapely/manual.html#prepared-geometry-operations
https://docs.djangoproject.com/en/dev/ref/contrib/gis/geos/#prepared-geometries


texts:
I'm trying to customize feature form with ui file, I followed https://nathanw.net/2011/09/05/qgis-tips-custom-feature-forms-with-python-logic/1 but the attributes values still null when I added feature.
This is ui file

So are there something I missed?


texts:
I have a shapefile with shapelength and shapearea, and would like to convert to WGS84 to use with ggmap.

Eventually, I would like to get plot google map of individual regions of the shapefile, and overlay it with other date

What can I use in R program to do this, i.e. convert to WGS84?


texts:
New to the site, and to geoJSON, and to python.

My goal is to read in a file with geoJSON data, and use it to generate shapes to be used in an application that is using python scripting.  Python has the ability to read it in, and do some manipulating.  I have found some examples, but they are not traversing the geoJSON hierarchy correctly.

I'm looking for a library, or instructions, on how to extract the geometry data from the geoJSON data file.  In other words, I would like to be able to say "get the list of features..." - then, for each of the features, I'd like to get the geometry, read the type, and the get the coordinates for that particular geometry item.  Repeat for each of the items in the list of features.

Any suggestions?


texts:
I can't get the syntax correct for a fairly simple if statement.  Right now it looks like this:

dim Output as Integer
If [HECTARES] >= 30 AND [HEIGHT] >= 17 AND [CROWN] >= 59
 THEN Output = 1
 ELSE Output = 0
END IF

I have no idea whether that's the correct use of "AND" but the expression wouldn't evaluate with only one expression (i.e. [HECTARES] > 30) either.

Thoughts?


texts: Failure to write LineString record to shapefile


texts:
I'm interested in converting a GeoDatabase to a MongoDB database.
more specifically, given I have a GeoDatabase (.gdb files), how do I convert these files to fit into a MongoDB Database (.bson files).

My tools are the ArcObjects SDK for .NET and MongoDB shell.

Thanks :)


texts: st_ delaunay returning empty geometry from 3D points

texts:
The Data

I have a 3D point dataset, one row looks like this

SELECT ST_AsEWKT (geom) FROM pc_processing.top_dist
SRID=32635;POINT(580992.318 4275265.0982 191.284)
... +1000 more

The problem

I want to triangulate in 3D so I try the

SELECT ST_AsEWKT (ST_DelaunayTriangles(geom,0,0))

But if I select the following I get

SELECT ST_AsEWKT (ST_DelaunayTriangles(geom,0,0))
FROM pc_processing.top_dist
SRID=32635;GEOMETRYCOLLECTION EMPTY
....

So, how do I return a Delaunay triangulation from 3D point data?


texts: Leaflet/HTML/GeoJSON: How to put a pop-up on a geoJSON point?

texts:
I don't have much knowledge of HTML/JS/GeoJSON, and have made a webmap from trial and error with many different tutorials. I have tried leaflet's tutorials but cannot get the desired results.

I have a functioning Leaflet webmap that shows GeoJSON points. I need to symbolize the points based on GeoJSON properties, and provide a pop-up when the point is clicked. I can use "pointToLayer" to turn the default blue markers into a custom point,

but when I try to incorporate a pop-up, the blue points turn back into default markers

How can I get a custom-symbolized point with a pop-up? What am I doing wrong?

Here is the HTML code I am using:

var geojsonMarkerOptions = {
radius: 5,
fillColor: "#33ccff",
color: "#000",
weight: 1,
opacity: 1,
fillOpacity: 0.9
};

$.getJSON("./geojson/kitchener_sidewalks2.geojson",
 function(data) {
    var geojson = L.geoJson(data, {onEachFeature: action_To_Perform_When_Marker_Is_Clicked_On_The_Map},
       {
        pointToLayer: function (feature, latlng) {
            return L.circleMarker(latlng, geojsonMarkerOptions);
        }
    });
    geojson.addTo(map);
    myGeoJSONLayers.addLayer(geojson);
 });

var overlayMaps = {"Kitchener Sidewalk Issues": myGeoJSONLayers};

L.control.layers(baseMaps, overlayMaps,myGeoJSONLayers).addTo(map);

function action_To_Perform_When_Marker_Is_Clicked_On_The_Map(feature, layer) {
if (feature.properties)
{
    var PopupText = [];
    PopupText.push("Accessibility Issue: " + feature.properties.Issue);
    PopupText.push("
Description: " + feature.properties.Descriptio);
    PopupText.push("
Recorded: " + feature.properties.CreationDa)
    layer.bindPopup("
" + PopupText.join("") + "
");
}
}


texts: SLD file with .png's to style points on a Geoserver OpenLayers map

texts:
The .png's are stored here: C:\Program Files (x86)\GeoServer 2.12.1\data_dir\styles.

I have seen many different ways to incorporate a .png into Geoserver and none have worked for me. I have no idea how to proceed.

These are a few formats that I have tried:

xlink:href="http://localhost:8080/geoserver/styles/Switch.png"
xlink: "Switch.png"
xlink:href="http://10.29.6.112:8080/geoserver/styles/Switch.png"

Does anyone have any advice or suggestions?

P.S. I am working from IP 10.29.6.112.


texts:
I want to doing geocoding with ArcGis geocoding service
but I receive this error:

How can I solve this problem ?


texts: Colour specific polygons on a postcode map - QGIS

texts:
I am working on a UK postcode map. It only has the first bit of the postcode, e.g. BR12. I have an excel spreadsheet with the longitude/latitude of specific points as well as their postcodes. Instead of plotting points I'd like to colour the postcode in which they are located.

Manually selecting the postcodes to colour would be a huge waste of time. Has anyone got a suggestion on how to proceed?

As per @Jospeh's suggestion, I have decided to go for the Rule-Based option on my .shp -however I don't understand how to add a rule as proposed below.


texts: Summing Specific Features of Overlapping Layers

texts:
I have several layers that are separated by categories (1-5 being 1 low importance and 5 high importance). Each layer only possesses 1 value.

I merged two layers simply using the Merge Vector function. Now I would like to create a new field which gives me the sum of categories. The final aim is to have an idea of which areas are more important (Those with the higher sum).

(Currently, the 0's in my merged layers are NULL values instead of 0's but I guess I can replace that with a simple xxx is NULL and then update selected features to 0)

**Desired Output**

| Layer A    | Layer B     |       Output |
|:-----------|------------:|:------------:|
|     0      |      2      |      2       |
|     1      |      0      |      1       |
|     1      |      2      |      3       |
|     1      |      0      |      1       |
|     0      |      2      |      2       |
|     0      |      2      |      2       |

How can I make this happen? Should I use Field Calculator?


texts: Keep Getting Black OSM Tiles in Openlayers

texts:
When trying to view OSM tiles using the external storage example of OpenLayers,
I keep getting black tiles when I zoom in and out.

I don't know what the issue is but every time I zoom in/out there are black tiles that appear.

Please help me.


texts: PostgreSQL select value where '100%' using bigint

texts:
I am trying to select data from a bigint column in PgAdmin using SQL:

SELECT FROM public.addressbaseplus WHERE uprn LIKE '100%'

However, I am getting the following error:

ERROR:  operator does not exist: bigint ~~ unknown
LINE 1: SELECT FROM public.addressbaseplus WHERE uprn LIKE '100%'
HINT:  No operator matches the given name and argument type(s). You might
need to add explicit type casts.
SQL state: 42883
Character: 47

Can someone explain to me what I am doing wrong and the error appears?


texts:
I am using in my project the Esri Ocean Basemap. I use the fullscreen control to enlarge it, but when enabled the map is not covering the whole screen.
The HTML structure is as follows:

The controls are loaded as
var map = new ol.Map({
    target: 'map',
    layers: [],
    controls: ol.control.defaults().extend([
        new ol.control.Zoom(),
        new ol.control.Attribution(),
        new ol.control.ZoomSlider(),
        new ol.control.MousePosition({
          projection: mercatorProjection,
            coordinateFormat: function (coordinate) {
                return 'Coordinates: ' + 'Lat ' + ol.coordinate.format(coordinate, '{y}', 3) + ' Long ' + ol.coordinate.format(coordinate, '{x}', 3);
            },
            target: 'coordinates'
        }),
        new ol.control.ScaleLine(),
        new ol.control.FullScreen({source:'full-size-container-custom'}),
        new app.Legend()
    ]),
    view: view

And the CSS is
.full-size-container-custom {
margin-bottom: 10px;
width: 100%;
height: 100%;
position: relative;
}
.full-size-container-custom:-moz-full-screen {
 height: 100%;
 background-color: #fff;
}
.full-size-container-custom:-webkit-full-screen {
  height: 100%;
 background-color: #fff;
}
.full-size-container-custom:-ms-fullscreen {
  height: 100%;
  background-color: #fff;
 }
.full-size-container-custom:-fullscreen {
 height: 100%;
 background-color: #fff;
}
map-container {
height:100%;
width: 80%;
padding: 0 0 0 0;
float:left;
}
.sidebar-map {
width: 19%;
height:100%;
padding: 0 0 0 1%;
float:left;
}
.map{
width: 100%;
height: 100%;
}
.map {
position: relative;
}

The aspect of the map can be seen in the following two screenshots. Do you have any suggestion?
Map in normal view

Map in fulscreen


texts: Finding polygons without right angles using Open Source GIS or ArcGIS for Desktop?

texts:
Now we are digitizing some buildings in specified area.

The obligatory rule for this work - in most cases buildings should have right angles.

We are using QGIS with CAD tools for this work but sometimes we make mistakes and create polygons with an irregular shape.

Does anybody know how can we find such polygons without right angles using open source GIS or ArcGIS?


texts:
I tried to change the projection of AMSR-2 in Python using:

'gdalwarp -geoloc -t_srs EPSG:4326 C:\GW1AM2_201301311114_050A_L1SGBTBR_1110110.vrt C:\test100.tif'

It said there are something wrong with 'EPSG:4326'(Invalid syntax). I do not know why?


texts: Autonumber field that is populate when a point is created

texts: I have a feature class called Ref_Point I need a field that creates a sequential number as points are added. Im using the feature class as a numbering tool for different jobs, so when I finish a job and clear the data I want the field to start back at number 1


texts: ArcGIS Server config to allow ArcGIS JavaScript printing?

texts: I have ArcGIS Server 10.21 on Windows Server 2012 R2, and have created a web application using ESRI JavaScript API v.3.16 that uses the dijit/Print that invokes the Export Web Map task on the server.  I have the Printing Tools service running on the server and can print any of the ESRI background maps just fine.  However, when I add any of my local feature layers, I get an error that it is unable to complete the process.  The error message in the PrintingTools log says that it "failed to create layer from service at http://myServer.myDomain.com/arcgis/rest/services/SampleWorldCities/MapServer/0".  If I replace the server URL with the local IP address (ex. 192.168.2.30) then it works fine.  Does anyone know what I am doing wrong here?


texts: Resetting all layers one time, when click is fired

texts:
Within a click function, I have the below code:

home.addEventListener('click', () => {

       myAPP.mapview.map.layers.items.forEach(item => {
            item.visible = false;
        });
       //........... }

However, after the button is clicked once, it seems to continuously run/execute. How could I assure it executes, i.e. 'resets' just one time, just when the button or click function occurs?


texts: st_astext generate an invalid wkt in postgis


texts:
I want to summarize fields per Group with python and need help to get unstuck. I'm trying to get the code below to produce a Dictionary of Groups and their total shape areas, like {'Group1':123456, 'Group2':654321}. The code below will produce a Dictionary with all Groups but with zeros for shape areas for all but the first Group who will get the shape area of only the first row in the feature class fcunion.

import arcpy
fcunion=r'D:\Folder\Database.gdb\fcunion'

fidfields = [f.name for f in arcpy.ListFields(fcunion, 'FID*')]
fidfields.extend(['field1', 'SHAPE@AREA','field2'])

#Create list of the groups
groups=[]
scursor = arcpy.da.SearchCursor(fcunion, fidfields)
for row in scursor:
    if row[5] not in groups:
        groups.append(row[5])
del row

#Summarize shape area for each group and create a dictionary
out_dict={}
scursor.reset()
for row in scursor:
    for g in groups:
        lst=[]
        if row[5]==g:
            lst.append(row[4])
        out_dict[g] = sum(lst)
del row, scursor

print out_dict


texts: Extending Polyline in ArcGIS without intersecting line?

texts:
I am trying to extend the ends of each line I have in a seismic survey 2 miles. I have tried using the extend line (Editor) tool but it needs an intersecting line to move to.

Does anyone know how to extend a previously created line a specified distance on both ends without intersecting another line?


texts:
So I have a shapefile with 160k elements.
It is in UTM and I calculated the area (should be meters, right?).
What I need now is the average size of each of my polygons.
Meaning, summing up a column and then divide this number by 160k.

Is there a straight forward way?

If there is no straight forward way, I'd rather do in Python than copying my attribute table to Excel and do it there.


texts:
I've been developing a Plugin on QGIS platform using PyQGIS. So far plugin does the following: You can select any feature of the already loaded vector layer in QGIS. After selection, plugin displays attribute names, and associated attribute values, explicitly for the selected feature. You can change those values, and save them.
Now I would like to add a possibility to define data type for every attribute that loaded layer contains. For example, I would like that first attribute named, let's say, ''road_id'' can only be 4-digit long integer, so if you enter string value or even a 5-digit integer, it gets rejected.
Is it possible to do that? If it is, which are the classes/methods I should be looking into?


texts: Prevent XML External Entity Injections in Geotools WMS Client

texts:
I'm using Geotools 11.5 and would like to protect against XML External Entity Injections when using the Geotools WMS client (gt-wms). The linked OWASP page suggests disallowing doctypes as a means of achieving this with a SAXParserFactory, which I believe is used by gt-wms. Is there any way of disallowing doctypes when using gt-wms?

Possibly related is that even in the later versions of org.geotools.xml.DocumentFactory:152-153 appear to have the lines of code needed to fix the issue, but they are commented out? The github history doesn't go back that far, can anyone shed any light on this?


texts: Fish passage connectivity analysis using geometric network traces?

texts:
We are looking at the connectivity between lakes and the ability of fish to pass between them, considering that downstream presents no barriers for these fish but swimming upstream does (barriers = steep gradients embedded in the stream network).

The best solution I have found is to build a geometric network, and run "trace upstream" and "trace downstream". My problem solving drops when I have to figure out how automate it for all 391 lake points, so that the output looks the same as you would get from running an OD Cost Matrix except instead of resistance values I would simply have "connected" or "not connected" or; even a table that resembles a spatial join where each unique waterbody ID has a corresponding record for each lake that it is connected to, as a fish is able to swim.

I have attempted using ModelBuilder and inserting iterate functions to loop through each record in the lake points feature class before the trace operation, but I am not getting the results that I would like (ie either it outputs the exact same feature class or it only outputs the last point). I am not highly experienced in ModelBuilder either so my sequence of functions may not be correct.

Edit Here is my model:


texts:
Is there a way to have a set of saved queries in QGIS?

I want to be able to bring in the same dataset, but with different values on different days, and be able to run certain queries on it with just a quick click rather than going through the query builder.


texts: Highlight the filtered Features using CQL_FILTER on map

texts:
I have a WMS service with multiple layers which I am consuming within a single source through Layers Param.

I need to filter features for few layers from the service. I am able to get the updated image with the filters from GeoServer , But the filtered features from the layers are not highlighted and it's difficult to locate the change on a dense map.

Below is my code to filter the features :

var multipleSrc = new ol.source.TileWMS({
    url: url,
    params: {
        'FORMAT': 'image/png',
        'LAYERS': getCommaSeperatedLayersFromHT(), //Ex : "Layer1;Layer2;Layer3;Layer4;Layer5;Layer6;Layer7;Layer8"
        'TRANSPARENT': 'TRUE',
        'CQL_FILTER': getCommaSeperatedLayersFromCQL() //Ex : "INCLUDE;INCLUDE;INCLUDE;INCLUDE;INCLUDE;INCLUDE;INCLUDE;Key='Value'"
    },
    serverType: 'geoserver',
    transition: 0,
    transparent: true,
    opacity: 0.5,
});

If I update the CQL_FILTER then the map is again regenerated with the revised filtered criteria.

Is there a way so that I can highlight the filtered features or I can indicate the filtered features in any way so that it's easier to identify them on map.


texts: MapProxy with WMTS source and reprojection

texts:
I am trying to consume Lantmateriet's WMTS source with EPSG:3006 and convert that to iOS TileOverlay compatible (GoogleMaps compatible) tile server. I cannot get it working. Can someone check this out and say what I am doing wrong.

mapproxy.yaml:

services:
  tms:
    use_grid_names: false
    origin: 'nw'
layers:
  - name: sweden
    title: sweden
    sources: [cache_sweden]
    tile_sources: [cache_sweden]

caches:
  cache_sweden:
    sources: [cache_sweden_original]
    grids: [GLOBAL_MERCATOR]
    meta_size: [4, 4]
    cache:
      type: mbtiles
  cache_sweden_original:
    sources: [source_sweden]
    grids: [grid_sweden]
    disable_storage: true
sources:
  source_sweden:
    type: tile
    grid: grid_sweden
    url: 'https://api.lantmateriet.se/open/topowebb-ccby/v1/wmts/token/mytoken/?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=topowebb&STYLE=default&TILEMATRIXSET=3006&TILEMATRIX=%(z)s&TILEROW=%(y)s&TILECOL=%(x)s&FORMAT=image/png'
grids:
  grid_sweden:
    srs: 'EPSG:3006'
    bbox: [-1200000.000000, 4305696.000000, 2994304.000000, 8500000.000000]
    origin: nw

This is what it looks now:

Capabilities of the WMTS source I am trying to use can be seen:
GetCapabilities.xml


texts: Creating and Analysing sales potential Territories in Carto

texts:
Is there any way to analyse territories and edit territories in Carto?

I do not need an algorithm to do it. I would need to aggregate polygons manually in order to create territories, see the results and keep tuning them until it makes sense.

I have polygons that represent neighbourhoods. Each polygon has data. Numbers of habitants, indexes, Number of companies, etc. I want to BE able to make/join polygons lets say to have 5 big areas. Each área should be identical so I must BE able to remove and add polygons to each area in order to understand if they are approximately similar.


texts: How to automatically sum up values within buffers?

texts:
I want to establish the population figures within multi-ring buffers around certain site. I have set up the buffers and am running a select by location query on each concentric ring.

Once the selection is made, I am manually accessing the table and using the Statistics function to find the sum of the population.

As I have 76 sites and 10 buffers per site, this is obviously taking a very long time to type up 760 individual results. Is there a way to automate this process and produce a table of the results?

I know nothing about Python scripting at the moment, but I would be interested to know if this could perform such an operation.


texts: Exporting selection to Excel file using ArcMap

texts:
I selected a list of facilities that intersect a polygon using ArcMap, and I want to export them to Excel.

Is there a way to do this without having to export the selection to a shapefile first?


texts: How can I display the symbology of the child layer of a parent layer?

texts:
I am working with layers in relations by a sub form (relations from 1 to n, relations established in the project properties named as "DB").
Parent layer: peoples with name (main form). Point layer (geopackage).
Child layer: the details of peoples fortune (sub-form without geometry geopackage).
How can I display the symbology of the child layer from the parent layer?
For example to display the totality of the fortune by proportional circles or an element of the fortune in particular like money cash (argent), buildings (immeubles), goods (marchandises).


texts: QGIS Network Analysis: QgsLineVectorLayerDirector.makeGraph(), attribute type QVector not available

texts:
relates to QGIS 2.14.6LTR on Win10, PyQt Version 4.10.2

When trying to makeGraph() on QgsLineVectorLayerDirector as described in the "QGIS Network analysis library" documentation (http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/network_analysis.html) like

director = QgsLineVectorLayerDirector(vl, -1, '', '', '', 3)
#[...]
psource = f_source.geometry()
ptarget = f_target.geometry()
#each assumed to be of type QgsPoint
#director has added a QgsDistanceArcProperter
#builder is a QgsGraphBuiler
#[...]
director.makeGraph(builder, [psource, ptarget])

I get an error message QgsLineVectorLayerDirector.makeGraph(QgsGraphBuilderInterface, list-of-QgsPoint): argument 2 has unexpected type 'list' See log for more details

Isn't the list internaly converted to QVector which is expected by makeGraph() https://qgis.org/api/classQgsLineVectorLayerDirector.html#af14ff34e31a2860edd156090c0f746ef

What should I pass to makeGraph() instead of my list?

ADD/UPDATE:
I tried

[...]
tp = QVector()
tp.append(psource)
tp.append(ptarget)
tiedPoints = director.makeGraph(builder, tp)
[...]

but name 'QVector' is not defined See log for more details. Isn't this imported with from PyQt4.QtCore import * ?

ADD: Obviously not:

Refer to also: https://stackoverflow.com/questions/39746241/pyqt4-convert-list-to-qvector


texts: ArcGIS error: Unable to save edits. Raster dataset already exists

texts:
When I edit my attribute table then try to save edits an error message appears on my screen: "Unable to save edits. Raster dataset already exist." Why does this error message appear?

My table of contents contains one feature class named "Buildings". Buildings have a lot of fields - building name (text), floor number (integer), picture of building (raster), etc. and their capacity about 100gb.

Until Thursday I have edited the attribute table and saved correctly but suddenly I received the above error message on Sunday. I edited some polygons without raster and saved without problem. However I try to edit polygons which contain a raster and then I want to save edits, but the error message does not let me.


texts:
I have the following VB code to produce a percentage in a field calculator.

" [FIELD_1] / [Total]*100 "

It works in VB. I've tried the following without luck.

!FIELD_1! / !Total! * 100

Can anyone translate this VB code to python?


texts: Calling Python script tool from ArcGIS model using relative path?

texts:
I am new to ArcGIS and facing a problem in calling a script tool from a model.

I have the following file structure:

E:\MyFolder

MyToolbox.tbx

ScriptTool
Model

Scripts

myscript.py

The ScriptTool uses 'E:\Scripts\myscript.py' as the script file and is configured to use relative path. Hence, even if I move MyFolder from E:\ to D:\ , the ScriptTool is able to find myscript.py and it executes successfully.

Also the ScriptTool is added to the Model. The Model works fine if MyFolder is in E:\ but if I move MyFolder to D:\ , the Model cannot run as it is not able to find the ScriptTool.

Please help me by suggesting how I can add the ScriptTool to the Model such that the model uses relative paths to the script tool.


texts:
I'm trying to load some XYZ data into a program called OpendTect. However, the orientation of the existing data is different from the grid orientation of the XYZ file. this makes it imposible to load the file correctly.

Is it possible to reorientate / rotate the XYZ grid using QGIS.


texts:
I have several .mbtiles with no (or little) overlapping areas and I would like to merge them into one file in order to be able to serve them more easily.

I'm looking for a command line tool that can do this.

tile-join seems to work only on vector tiles. When I run to on my raster tiles it crashes:

tile-join first.mbtiles second .mbtiles -o combined.mbtiles
libc++abi.dylib: libc++abi.dylib: libc++abi.dylib: libc++abi.dylib: libc++abi.dylib: terminating with uncaught exception of type protozero::unknown_pbf_wire_type_exception: unknown pbf field type exceptionterminating with uncaught exception of type protozero::unknown_pbf_wire_type_exception: unknown pbf field type exceptionterminating with uncaught exception of type protozero::unknown_pbf_wire_type_exception: unknown pbf field type exceptionlibc++abi.dylib: terminating with uncaught exception of type protozero::unknown_pbf_wire_type_exception: unknown pbf field type exceptionlibc++abi.dylib: terminating with uncaught exception of type protozero::unknown_pbf_wire_type_exception: unknown pbf field type exception
Abort trap: 6

gdal_merge.py works on .tiff files but I have .mbtiles as input.

Update:

I managed to convert the .mbtiles to .tiff using

gdal_translate input.mbtiles output.tiff
gdal_translate input2.mbtiles output2.tiff

Then merge them using

gdal_merge.py -co COMPRESS=LZW  *.tiff

This gives me one out.tiff but then how to get a .mbtiles back? goal_translate gives me this error:

gdal_translate -of mbtiles out.tif  out.mbtiles
Input file size is 81920, 81920
0ERROR 6: GDALDriver::Create() ... no create method implemented for this format.

any ideas?


texts: Trigger AutoUpdaters in Stand Alone Python Script ArcFM

texts:
I am trying to insert records from an FGDB featureclass into an "ArcFM-ified" SDE(Oracle 11g) feature class with a stand alone python script.  I want the AutoUpdaters on that feature class to fire.  I have downloaded and looked at the examples on Schneider's exchange site.  I have also seen the few other questions on here like Objects in this class cannot be updated outside an edit session inside ArcFM?

According to what I have read I need to checkout an ArcFM license, set the

mmAutoUpdaterMode = mmAutoUpdaterMode.mmAUMStandAlone

and when I run the script the auto updaters should fire.

I run my script, the records are successfully inserted into the feature class, but no AutoUpdaters are fired.

Why are they not firing?

I am using arcpy.InsertCursor because I could not get arcpy.da.InsertCursor to work and it was suggested in this post on Schneider's exchange site.

Here is my code:

import arcpy
import os
from arcgisscripting import create
from win32com.client import Dispatch
gp = create(10.1)

app = Dispatch('Miner.Framework.Dispatch.MMAppInitializeDispatch')
au = Dispatch('Miner.Framework.Dispatch.MMAutoupdaterDispatch')
runtime = Dispatch('Miner.Framework.Dispatch.MMRuntimeEnvironmentDispatch')

from enumerations import mmRuntimeMode,mmLicensedProductCode,mmAutoUpdaterMode

runtime.RuntimeMode = mmRuntimeMode.mmRuntimeModeArcServer
app.Initialize(mmLicensedProductCode.mmLPArcFM)
au.AutoUpdaterMode = mmAutoUpdaterMode.mmAUMStandAlone

sde_conn = r"C:\Users\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\GDB.sde"
target_fc = os.path.join(sde_conn, 'Schema.FeatureDataset', 'FeatureClass')
source_fc = r'C:\FGDB.gdb\FeatureClass'

feature_list = []
fields = [
          'RECORDID',
          'SHAPE@',
          'SUBTYPECD',
         ]

def populateList():
    with arcpy.da.SearchCursor(source_fc, fields) as scur:
        for row in scur:
            feature_list.append(row)

def sampleArcFM():
    # code for geoprocessing goes here
    editor = arcpy.da.Editor(sde_conn)
    editor.startEditing(False, True)
    editor.startOperation()

    icur = arcpy.InsertCursor(target_fc)
    for rec in gen_list:
        row = icur.newRow()
        rec_id = rec[0]
        shape = rec[1]
        subtype = rec[2]
        row.setValue("RECORDID", rec_id)
        row.setValue("SUBTYPECD", subtype)
        row.shape = shape.firstPoint
        icur.insertRow(row)

    editor.stopOperation()
    editor.stopEditing(True)

    del icur

    app.Shutdown
    del app, runtime, au

populateList()
sampleArcFM()


texts:
Is it possible to digitize a polygon on qgis in freehand?
I would like to create a polygon vector without clicking on each point with the mouse.


texts: Python to create new fields in many layers / QGIS

texts:
I have a number of line layers on which I need to perform the same task: add 4 columns with the start and end points of each line.

The expressions that I use to do this are:

$xat(0)
$xat(-1)
$yat(0)
$yat(-1)

However, I must create a new column and update it for every layer, but I am sure that it can be automated.   Any ideas on how to perform the creation of 4 new columns on every layer in the layers panel in QGIS?


texts: ArcGIS 10 (SP 3) - Geometric Network - which features are moved?

texts:
In ArcGIS 9.3 Geometric Network, you can select which features can be moved when snapping.

In ArcGIS 10 (SP 3) this functionality has disappeared and you either choose to have snap tolerance or not.
Does anyone know why this functionality is no longer available? When snapping 3 feature classes to create a geometric network, how do we know which feature class has moved? We only want 1 of the 3 feature classes to move.
Does anyone know how we can define which feature class is to move (and snap) to the other feature classes?