texts:
I am running a script to export a series of mxds to pdf via ArcPy in an ArcToolbox tool (code below).
I already reduced it to the minimum of code necessary but still, NON REPRODUCIBLE, I get a Microsft Visual C++ Runtime Error r6205 "pure virtual function called". This happens with different mxd on all kinds of different Win7 machines. ArcGIS version is 10.1.

Here is the code:

#Load necessary modules
import arcpy, os, string

#Read input parameters from script tool
mxdList = string.split(arcpy.GetParameterAsText(0), ";")
outPDFpath = arcpy.GetParameterAsText(1)
resol = arcpy.GetParameter(2)
imQal = arcpy.GetParameterAsText(3)
comVec = arcpy.GetParameterAsText(4)
comImg = arcpy.GetParameterAsText(5)
picSym = arcpy.GetParameterAsText(6)
conMark = arcpy.GetParameterAsText(7)
embFo = arcpy.GetParameterAsText(8)
comJPG = arcpy.GetParameterAsText(9)

mxdList = [mxdPath[1:-1] for mxdPath in mxdList]

#Statusbar for rough prediction of duration
#totLength = len(mxdList)
#arcpy.SetProgressor("step", "Exporting Maps to PDF...",
#                    0, totLength, 1)

#If folder not existing -> create
if not os.path.exists(outPDFpath):
    os.makedirs(outPDFpath)

# Send Message to Operator
#if totLength <= 10:
#    arcpy.AddWarning("\n Dear Operator, this might take a while. Get an icecram or a quick coffee. \n")
#if totLength > 10:
#    arcpy.AddWarning("\n Dear Operator, this might take a damn good while. Get a beer or two. \n")

#Loop through each MXD in the list, export, output PDF
for mxdPath in mxdList:
    mxd = arcpy.mapping.MapDocument(mxdPath)
    mxdName = os.path.splitext(os.path.basename(mxdPath))[0]
    PDFPath = outPDFpath + "\\" + mxdName + ".pdf"
    #Short message to the operator which Map is being processed now
    #arcpy.SetProgressorLabel("Exporting Map to PDF:\n" + mxdPath)
    arcpy.mapping.ExportToPDF(mxd, PDFPath,
                              resolution = resol,
                              image_quality = imQal,
                              compress_vectors = comVec,
                              image_compression = comImg,
                              picture_symbol = picSym,
                              convert_markers = conMark,
                              embed_fonts = embFo,
                              jpeg_compression_quality = comJPG)

    # Update the progressor position
    #arcpy.SetProgressorPosition()
    # Send Message to Statuswindow
    #arcpy.AddMessage("\n Supersuccess: " + PDFPath + " created! Processing next file. \n")

#Remove variable reference to file
del mxd, mxdList, outPDFpath, mxdName, PDFPath

Previous research indicated that all status message related functions might be likely to cause such a crash, so they are commented out.

Does anyone see an error in my code? I am just stuck with finding a solution.


texts: Writing value into specific column in ArcGIS table?

texts:
I am writing a Python script to add fields and write data into the ArcGIS table. I can successfully edit and delete columns or fields in the table. But I am stuck at writing data into the fields. I got a code by which we can get the type and others info about the fields. But I was unable to write something on that particular field.

Here is the script I got from their forum -

import arcpy

feature_class = "c:/data/counties.shp"

# Create a list of fields using the ListFields function
fields = arcpy.ListFields(feature_class)

# Iterate through the list of fields
for field in fields:
    # Print field properties
    print("Field:       {0}".format(field.name))
    print("Alias:       {0}".format(field.aliasName))
    print("Type:        {0}".format(field.type))
    print("Is Editable: {0}".format(field.editable))
    print("Required:    {0}".format(field.required))
    print("Scale:       {0}".format(field.scale))
    print("Precision:   {0}".format(field.precision))

The picture below is the table I have

My Primary Key is COPROPCD and I want to change the other fields based on my other dataset which has the common primary key.

Update 2
I tried the following code but it still doesn't work and get errors.

import arcpy
 feature_class = "c:/BaseMapsLocalCopy/Export_Output.shp"

 with arcpy.da.UpdateCursor(feature_class, ['COPROPCD','LEGAL']) as cursor:
     for row in cursor:
         #Check if "COPROPCD" exist as a key in the dict
         #if row[0] in dataset.keys():
         #logic for updating the dable
         #This will upload the field2 column to the value of the key
         if row[0] == 'R040222':
             row[1] = 'GORU'
             cursor.updateRow(row)

And the error it is giving is that there no column that was specified. The exact error is as below -

Final Update
It works now. My file was wrong. Correcting the feature_class solved my issue.


texts: Shifting display of overlapping lines in QGIS?

texts:
When points are overlapping, there is this property which allow to automatically display the lot of them separately around where they are, called 'Point displacement'. But it doesn't work for lines, even so it seems to me quite conceptually feasable in order to achieve something like that :

I absolutely need to see the differents lines which in reality are all in the same place (I am working in telecomunication networking). The only way I see for now is to really create different lines like in the picture above, thus creating spatial mistakes.

I'm using QGIS 2.14.


texts: Extract the north-west point of a shapely rectangle

texts:
I am using Python in a Jupyter notebook.
I'have created a geometry with shapely and I need to extract the coordinate of the north-west corner of it.
from shapely.geometry import Polygon

lng = -180
lat=-90
lon_point_list = [lng, lng+10, lng+10, lng]
lat_point_list = [lat, lat, lat+10, lat+10]

polygon_geom = Polygon(zip(lon_point_list, lat_point_list))

of course this shape is an example but my shape is always a rectangle
Is there a trick to perform it ? (I can use another library such as geopandas)


texts: Calculating distances (differences?) between two regions?

texts:
I am building a machine learning system which predicts the location of tweets as a bounding box (using the google maps API among other things). Since the twitter API also supplies a bounding box, I am trying to calculate the distance between my predicted box and the actual one.

Up to now, I have just been calculating the center of each box and using the distance between them as my measure, but this method is problematic because the regions can overlap partially or be contained within each other, etc. Is there a better way to determine the difference between the two regions?

The only idea I've had so far is to measure the overlap as a percent of each box's area and then maybe use that as a weight somehow, but GIS is not my field, so I know very little about any existing techniques.

My other idea was to model this as more of a geometry problem by just looking at this as the average distance between two random points selected respectively from each region. This obviously doesn't take into account the population differences between (and inside of) the regions, but for small regions this would be valid. This method would eliminate the problem of potentially overlapping regions.


texts: Can't get pyshp to work when converting csv to point file

texts:
I've been hammering at it for ages, but I just can't get pyshp to work. The CSV file looks like this:

lon,lat
-122.5163113, 45.55324148
-95.1719409, 29.97979461
153.0072888, -27.45988559
...

There are many more rows to be had. Using a lot of examples on the internet I've written the below code, but whenever I try to run it I get the error:

Traceback (most recent call last):
  File "H:/SU5050/shapefile.py", line 1, in <module>
    import shapefile
  File "H:/SU5050\shapefile.py", line 14, in <module>
    w=shapefile.writer(shapefile.POINT)
AttributeError: 'module' object has no attribute 'writer'

I've uninstalled and reinstalled pyshp (using pip) and when I try to simply import shapefile it doesn't error out, just when I try to run it.

import shapefile
import csv
#LON IS X (0)
#LAT IS Y (1)
file_location = 'H:\SU5050\WORK.csv'
out_file = 'H:\SU5050\GPS_Pts.shp'
x,y=[],[]
with open(file_location,'rb') as csvfile:
    r=csv.reader(csvfile,delimiter=',')
    for i,row in enumerate(r):
        if i>0:
            x.append(float(row[0]))
            y.append(float(row[1]))
w=shapefile.Writer(shapefile.POINT)
w.field('x','f',10,8)
w.field('y','f',10,8)
for j,k in enumerate(x):
    w.point(k,y[j])
    w.record(k,y[j])
    w.save(out_file)

Anyone have any idea what I'm doing wrong?


texts: ERDAS Imagine batch wizard edit replacement pattern help

texts:
I'm using ERDAS Imagine to batch produce NDVIs from NAIP imagery.  I need to also replace the default file name (e.g. "m_3610132_nw_14_1_20100503") with a user friendly file name (e.g. 3610132nw).  How can I change the edit replacement pattern in the batch wizard to output a subset of the original filename?

Example:
z:/project/studyarea_8_naip/ok_naip/temp_out/$(Input.pathright)$(Input.root)

yields => m_3610132_nw_14_1_20100503

However, I want the output to be simply 3610132nw


texts: Subsetting mutliband raster in QGIS using OTB BandMathX

texts:
I am attempting to do a simple extraction of a multiband (360 bands) raster in QGIS. I have come to the conclusion I need to use BandMathX from OTB, since I am interested in generating a new multiband raster with the values from the previous, above a certain value, in this case -1.5.
I have tried using the following in QGIS (this is copypasted from the log-file):
{ 'il' : ['dir/file.tif'], 'out' : 'TEMPORARY_OUTPUT', 'exp' : 'im1>{-1.5}', 'incontext' : '', 'outcontext' : 'TEMPORARY_OUTPUT', 'outputpixeltype' : 5 }

2021-03-31 11:18:48 (INFO) BandMathX: Using expression: im1>{-1.5}

2021-03-31 11:18:48 (FATAL) BandMathX: itk::ERROR:

ParserXImpl(0000021924F14460):

Message: Can't evaluate function/operator ">": Internal error.

Formula: im1>{-1.5}

Token: >

Position: 3

I have read the section on BandMathX in the OTB CookBook and I cannot see how the expression is in valid.
I have also tried using the same expression without the '{}'.
I hope some of you guys can help me out. I have already been searchin for similar posts, but were unable to find any.


texts: Avoiding lost layers when using QConsolidate?

texts:
I want to save all the files in my project and the project file locally on my computer. The files are currently saved on a external hard drive.

I have tried QConsolidate before, and it worked, but this time I get a lot of lost layers when trying to reopen my project from my local disk?

I am using QGIS 2.12 Lyon and my files are *.shp and *.tif.


texts: Calculate polygon intersection without boundary intersections

texts:
I'm trying to intersect the shapefiles for US congressional districts and US counties. I can intersect manually (Vector->Geoprocessing Tools->Intersect) or with python (processing.runalg('qgis:intersection', input, input2, output)). Here's a picture of my goal:

My problem is that my output file seems to have rows for each boundary intersection, eg. between the left triangle and the square. I get duplicates in the output file, where the area of the intersection is really tiny, and sometimes a county has the state name of its neighboring state (with which it shares a border).

I don't want to keep line intersections like this, only the polygon intersections. How can I do this in QGIS?


texts: Downloading Scotland's national power grid from OpenStreetMap

texts:
I've been tasked with finding a complete dataset of Scotland's electricity grid network. This comprises a series of cables carrying 132 to 400 kV. I'm assuming that many of these cables are above ground. The National Grid publish an open dataset of both cables and pylons for England and Wales. However the grid in Scotland is managed by Scottish and Southern Energy and Scottish Power, and neither publishes GIS data.

I've noticed that cables and pylons are shown on OpenStreetMap, but only draw when you zoom in beyond 1:10,000 or so. I don't understand the complexities of OSM enough to make the most of the tagging system, so I could well be missing something fundamentally important here. I downloaded the 'Scotland latest' dataset from Geofabrik but this seems to be generalised, and doesn't include the cables I need.

I have had moderate success in QGIS importing lines for a trial area using the Vector -> OpenStreetMap -> Download Data route. However I cannot download national coverage this way.

I've also installed Java OpenStreetMap Editor, on the back of information online that suggests you can use this to download data. Unfortunately I've not found a way to query data for download using tags and you cannot download national-scale datasets.

Rather than take a 'dump' of all OSM layers, I would quite like to be able to state the layers I require and, say, a bounding box for which I wish the data within.

Ultimately the grid would be used in ESRI ArcMap in a map for publication.


texts: Longest common prefix of geohashed points

texts:
I currently use Geohash to index some points and retrieve information faster by providing a geohash prefix for a given point for a given distance..

For example if I get a given point at '1600 Pennsylvania Ave, Washington, DC' with geohash dqcjqcpsr4t0 and I get all point that geohash starts with dqcjqc.

I would like to know if there is an algorithm and index solution to get the points  ordered by the longest common segment without iterating through them all.

For example: I got a given point at '1600 Pennsylvania Ave, Washington, DC' with geohash dqcjqcpsr4t0 and I would like to get the nearest point first:

point with geohash dqcjqfzg038b will be returned first (common prefix 'dqcjq'=>5 letters)
point with geohash dqcjz4p344cu will be returned next (common prefix 'dqcj'=>4 letters)

The idea behind this is to get the nearest points first by providing the full (8 chars) geohash of a given point (in this case we don't use geohash to get points for a given box)


texts: Postgis ST_MAKELINE change dimension of points

texts:
I have multiple 2D points(lng=X lat=Y altitude=Z==NULL timestamp=M) mixed with multiple 3D points (lng=X lat=Y altitude=Z timestamp=M). I am going to apply ST_MAKELINE (Postgis) to those points. Then I realized that function auto fill 0 to Z axis of all 2D points, but I expected NULLs. So is there any function to avoid this?

Probably my description was incorrect. I am trying to compress multiple mixed dimension points with ST_SIMPLIFY. before I call ST_SIMPLIFY I need to make all points into a LINE via ST_MAKELINE. Then ST_MAKELINE auto convert 2D points into 3D point which is not expected. NULL should be NULL rather than 0.

Pseudo query: ST_SIMPLIFY(ST_MAKELINE([set of points]))

Update: What I have done so far (awaiting better solution)

I made a copy of mixed points as raw_points and then convert all mixed points into 2D (basically drop Z) then call ST_MAKELINE with all 2D points. Then call ST_SIMPLIFY to get a set of compressed points as compressed_points

Finally recover altitude information by following

Pseudo code:

select raw_point
from compressed_point
 LEFT JOIN raw_points on ST_M(compressed_point) = ST_M(raw_points);


texts: Exporting output in ASCII from a QGIS model

texts: I'm trying to find a geoalgorithm in QGIS's Graphical Modeler that can convert and export a vector to ASCII. Could anyone enlighten me if there's an option for that? Perhaps through some script but I really don't know where to start. I'm guessing one that uses some translation tool from GDAL.


texts: What are good sources to find GIS scholarships?


texts: Accuracy issues between data in QGIS and Google Earth

texts:
Wondering if anyone can answer this question. I have a kml in Google Earth. When I import it to QGIS, and then import the Google Earth images they are in slightly different places. I used a QGIS plugin so maybe that is the issue. Potentially the plugin imports the terrain images incorrectly.

Google Earth - Original
Zone 18M
209551.00 m E
9131622.00 m S

QGIS - First import (note, Zone 18M is not an option in QGIS...potential answer?)
209581.516575, -868393.873839
*I think the incorrect Y coordinate has to do with zone datum points, which typically have large offsets if I remember right. But look at the X coordinate, its 30 meters off! This is the main part of my question. Why are the x and y coordinates changing?

QGIS - Coordinates given at bottom of screen
-77.633805064386,-7.847746027456, which when plotted on Google earth line up where the original google earth point was.

One last issue. I imported water data from divaGIS for rivers, and the river is definitely a km or so to the east of the google earth image river. Should I trust divaGIS data? They have a suspicious quantity available for free.

http://www.diva-gis.org/gdata

Thanks!


texts: Getting Annotation Text Geometry as Polygon Using ArcObjects?

texts:
How can i get the rectangle boundary of an annotation object using ArcObjects?

Once we select an object on the map, there is a boundary as shown in the first picture.

Also is there an option to get the minimum rectangle that covers the text as shown in second picture?

I tried

    IAnnotationFeature annoFeature;
    annoFeature.Annotation.Geometry;

but this only gives me a polyline geometry which seems to be the baseline of text


texts: Using boundaries of city to cover world outside of it in solid color?

texts:
I need to create a mask to cover the rest of the world outside of a single city. (Here are Boston boundaries to use as reference for the boundaries of a city).

I've been told by support to open those boundaries "create a new geospatial file that would be made up of a rectangle that covers the whole earth with a hole in the shape of your city in the middle."

From there, I need to upload that file to Mapbox Studio and add it to cover the rest of the world outside of my city. I'm looking to style OpenStreetMap in Mapbox afterwards for that city. So this is a workaround for not being able to edit a single city in Mapbox Studio.

How would I create this mask?


texts:
I recently came across this new map projection, the Authagraph. It seems pretty interesting and I'd like to play around with it, but I can't seem to find the math for projecting it, short of deriving it from scratch.

Is the math behind the Authagraph documented anywhere online?


texts: Error while importing shapefile using shp2pgsql: Coordinate values are out of range

texts:
I have a shapefile with 180 points. I'm using shp2pgsql-gui via pgAdminIII. This is an error it gives me when I try to import shapefile.

Failed in pgui_exec(): ERROR:  Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

This is information about coordinate system:

Geographic Coordinate System:   GCS_WGS_1984
Datum:  D_WGS_1984
Prime Meridian:     Greenwich
Angular Unit:   Degree

Can someone help me?


texts:
I am developing a print map function in an OpenLayers 3 application it includes some functionality such as if its polygon then area of that polygon will be display OR if its line length will be display.

At present only line or polygon is displaying but its popup like feature length or feature area its not displaying

Please check attached screenshots.

After clicking on print button only map get printed as follows

But I want print of this with area only that feature is displaying


texts: Extracting surface points in module lasground of QGIS?

texts:
I need an advice. I have Lastools installed in QGIS 2.18.16 (Las Palmas). But the problem is that I cannot extract LIDAR from the initially classified data, according to the instructions on the official site https://docs.qgis.org/2.8/en/docs/training_manual/forestry/basic_lidar.html The manual is like a different lasground module. Lastools, specifically the lasground module is completely different for input and output data. Check image.

Can you tell me how to enter the data in the lasground module to extract ground data?


texts: Python Postgres Multiprocessing speed increase

texts:
I'm trying to speed up my PostGIS queries using multiprocessing. My current setup is using python and psycopg2 such as below. Although this has given a speed increase, it still seems like there are bottlenecks preventing rapid speed increases and I'm not sure where to go next.

I've increased a lot of the postgres parameters as suggested in 'Performance Tuning Postgres', but when I run this in AWS I never seem to be anywhere near maxing out RAM or I/O which should be the limitations for DB activity apparently. Can anyone suggest other methods for speeding this up?

> import os, sys, psycopg2, multiprocessing, time
> start = time.time()
> conn = psycopg2.connect("dbname=template_postgis_20 user=postgres")
> cur = conn.cursor()
> cur.execute("""SELECT count(*) FROM sites""")
> count = cur.fetchall()
> def getOidRanges(rownums, count):
>
>     oidranges = []
>     for row in rownums:
>         minoid = int(row[0])
>     return oidranges
>
> def mp(rownums, whereclause):
>
>     for row in rownums:
>         if row[0] == whereclause:
>             gid1 = int(row[0])
>             cur.execute("""
>             UPDATE sites SET postcode=(SELECT field62 FROM (SELECT field62, COUNT(field62) FROM addressbaseplusbh1_2
>             WHERE ST_Within(addressbaseplusbh1_2.geom, (select geom from sites where gid={0})) GROUP BY field62 ORDER BY count DESC)
>             as postcode LIMIT 1) WHERE gid = {0};""".format(gid1))
>             conn.commit()
>
>     return
>
> if __name__ == "__main__":
>
>     nbrquery = ("""SELECT gid FROM sites ORDER BY gid;""")
>     cur.execute(nbrquery)
>     rownums=cur.fetchall()
>
>     cores = (multiprocessing.cpu_count()-1)
>     procfeaturelimit = 1
>
>     oidranges = getOidRanges(rownums, count)
>
>     if len(oidranges) > 0:
>         pool = multiprocessing.Pool(cores)
>
>         for oidrange in oidranges:
>
>             whereclause = oidrange[0]
>             jobs = pool.apply_async(mp, (rownums, whereclause))
>
>         pool.close()
>         pool.join()
>         jobs.get()
>
>         try:
>             conn.commit()
>             cur.close()
>             conn.close()
>             end = time.time()
>             print end - start
>         except:
>             pass

EDIT:

@Craig, would it work then just to have this as the executed block?

curs.execute("""UPDATE sites SET postcode=(SELECT field62 FROM (SELECT field62, COUNT(field62) FROM addressbaseplusbh1_2 WHERE ST_Within(addressbaseplusbh1_2.geom,
(select geom from sites where gid={0})) GROUP BY field62 ORDER BY count DESC)
as postcode LIMIT 1) WHERE gid = {0};""".format(whereclause))
return


texts:
I want listen to Edit Session in ArcMap. But every time when i click on my button in ArcMap, Applicaation crash down.  If i debug,  it stopped at AddHandler Events.OnStartEditing, AddressOf Events_OnStartEditing under Sub Button1.

Error:
Task    In order to evaluate an indexed property, the property must be qualified and the arguments must be explicitly supplied by the user. ESRI.ArcGIS.Editor.IEditTask

Code:

Imports System.Runtime.InteropServices
Imports ESRI.ArcGIS.esriSystem
Imports ESRI.ArcGIS.Editor
Imports ESRI.ArcGIS.Desktop.AddIns
Imports ESRI.ArcGIS.Geodatabase
Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.Framework
Imports ESRI.ArcGIS.ArcMap
Imports ESRI.ArcGIS.SystemUI

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

    Public Sub New()

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

    Protected Overrides Sub OnClick()
        AddHandler Events.OnStartEditing, AddressOf Events_OnStartEditing
        AddHandler Events.OnStopEditing, AddressOf Events_OnStopEditing
    End Sub

    Private ReadOnly Property MyEditor() As IEditor
        Get
            Dim m_editor As IEditor
            Dim editorUid As New UID()
            editorUid.Value = "esriEditor.Editor"
            m_editor = TryCast(My.ArcMap.Application.FindExtensionByCLSID(editorUid), IEditor)
            Return TryCast(m_editor, IEditor)
        End Get
    End Property

    Private ReadOnly Property Events() As IEditEvents_Event
        Get
            Dim m_editorEvents As IEditEvents_Event = TryCast(MyEditor, IEditEvents_Event)
            Return TryCast(m_editorEvents, IEditEvents_Event)
        End Get
    End Property

    Private Sub Events_OnStartEditing()
        MsgBox("you are now editing")
    End Sub

    Private Sub Events_OnStopEditing()
        MsgBox("You will stop editing")
    End Sub
End Class


texts: Geoprocessing tool- output file path location for personal geodatabase in ArcGIS Server

texts:
I want to know about using a geoprecessing tool:

How do I define an output file path location for a personal geodatabase in an ArcGIS Server geoprocessing folder path?


texts: Making Feature to Raster tool run faster?

texts:
Suppose you have a feature class that you need to convert to a very detailed raster (where each cell size is quite small). With a cell size of 0.0092908447 in the Feature to Raster tool, the procedure takes hours to complete (mine was running for 4-5 hours with some feature class never seeming to end).

How can one make this run faster?


texts: Converting geometry of elements of table to move it from one srid system to another

texts:
I have a table containing a column geom loaded from a shape file which was defined in lambert II system (27582), but when I loaded it I specified that I want to work in lambert 93 (2154).
Unfortunately when I load it in QGIS the coordinates are always in lambert II.
Is there a function allowing me to convert these coordinates and make them pass from lambert II (27582) to lambert 93 (2154)?
I would like to work only in lambert 93.


texts:
A similiar question has been asked recently as How to make GRASS and SAGA tools capable to use in QGIS？

However I have some specific queries and process I have tried which are probably better in a new question for clarity.

The process I used is as follows:

Installed QGIS using OSGeo4W. Installed Dufour 2.1 under Program files. Apps directory shows saga_gui.exe (2.0.8). Grass plugin working ok. Install looks ok.
The SAGA instructions mention a configuration dialogue without being a bit more accurate. https://www.qgis.org/en/docs/user_manual/processing/3rdParty.html. I assume this is in QGIS Settings \ Options \System variables, and that you should be able to put in the path to the .exe? If so I cannot seem to get this to work. Or maybe it's somewhere else.
Also mentioned is that SAGA 2.1 is needed. I've downloaded this and, as it appears to be a "standalone programme", put the whole thing into the Dufour folder rather than reinstalling everything with the aim of configuring the path to the .exe from QGIS. Is this correct, as usually installs put files all over the place?

From my quick look at the program SAGA seems to look pretty good. However there appears to be only basic file support (dxf, txt, etc) and no ability to directly load MapInfo files which is what I was interested in.

How do I run it from QGIS?


texts:
I am writing a program in python that converts raster files using gdal, and using VRT as an intermediate format.  One of the things that i need to do is convert bands between different formats.  I am pretty sure that you do this using Pixel Functions.  Since python doesn't have this functionality, i have to use C++ for this part.  I used some code from the VRT tutorial website, but it doesn't seem to work.

GDALAddDerivedBandPixelFunc("MyFirstFunction", TestFunction);

typedef CPLErr
(*GDALDerivedPixelFunc)(void **papoSources, int nSources, void *pData,
                        int nXSize, int nYSize,
                        GDALDataType eSrcType, GDALDataType eBufType,
                        int nPixelSpace, int nLineSpace);

I implement that function with the code from http://www.gdal.org/gdal_vrttut.html, but it doesn't seem to be working.  I am calling gdal_translate to do the conversion process, but it seems like the pixel function does not stay.  The addpixelfunc returns a result of 0.  Is there some way to write a GDAL extension to have this function be permanent?

Update:

I have been able to correctly create a pixel function and get it to work if i manually convert the file in the same C++ program as the GDALAddDerivedBandPixelFunc resides, however i am still interested in how i could make this permanent. Although it may not be needed since i might want to let the user change some numbers in the function.


texts: Is there a way to make multiple ring buffers into the raster format? The distances are 0.25,0.5,0.75 & 1. Since they are less than 1 it shows as an error while using the polygon to raster tool. Eventually, I want to reclassify them as ranks & eventually use weighted overlay for the analysis. (ArcGIS Pro)


texts: How do you change starting point on a polyline? (split tool)

texts: There must be an easy fix for this.  I am using the split tool and I want it to split a polyline into equal parts.  For some reason, it is starting in the middle of my polyline.  How do I force it to start where I want it to? I am using ArcGIS v. 10.1.


texts:
I have a set of raster area results, that represent a prediction for a specific loccation.  I would like to measure accuracy of two methods, given the actual values, estimate-A method values and estimate-B method values.

At the moment, for each raster pixel, I'm getting the difference abs(ActualLayer.pixel.value - EstimateLayerA.pixel.value), and using the average difference for a given area as a measure of accuracy (closer to 0 would be better).

This is ok, but are there better methods to quantize the accuracy of estimates against the actual for a given area?


texts: ArcPy Select by attribute using create date

texts:
I'm writting a Python script that will look for newly created features and update them. This script will be automated to run every day at the end of the day.
I'm using the Select by attribute tool to select features created on or before now (variable). I'm using UTC time zone because the "Create Date" field was created by the Enable Editor Tracking tool and by default the date is in UTC time zone.
import arcpy, datetime
from datetime import timezone

Now = datetime.datetime.now(timezone.utc)
arcpy.SelectLayerByAttribute_management(LayerName, "New_Selection", "CreateDate <= Now")

I'm getting an invalid expression error
arcgisscripting.ExecuteError: ERROR 000358: Invalid expression
Failed to execute (SelectLayerByAttribute).

How do I use the variable Now in the expression?


texts:
The problem I'm dealing with is the old 180 meridian problem, except I'm not trying to plot or render shapes that cross it, I need to do intersections of lines that cross it. The only spherical function available in geotools appears to be orthodromic distance, so I need to do the intersections in 2d.

Can I take the naive approach and simply copy the WGS83 WKT and change the PM to 180? I can't use the UTM zones because my lines are long (airport to airport, i.e. Hong Kong to Vancouver).


texts: How to fill the column in attribute table with same values in QGIS?

texts:
I have about 19000 records, I added new column and I would like to fill it with  string values (eg. brand)- same value for all rows.

I'm not considering copy/paste. I know I could do it with exporting to .csv or something and edit it excel, but that feels a bit... lame and unmethodical?
So, any smart solutions?


texts: Going from hillshade to percent shaded using gdaldem?

texts:
I assumed hillshades output percent shaded (scaled to 255) at the given azimuth and altitude but this is incorrect. A flat raster at 80 degrees altitude has a maximum value of 251 instead of 255. How can I go from a hillshade to percent shaded and can I simply divide by the cosine of the 10 degrees in this case to get the correct value or is it more complicated?

Here is some R code with gdaldem that may help you see what I'm seeing:

library(raster)
r <- raster(ncol=100,nrow=100)
values(r) <- rep(1,ncell(r))
writeRaster(r,"flatelevation.tif")
system("gdaldem hillshade -az 180 -alt 80 flatelevation.tif hillshade_flat.tif")
hs <- raster("hillshade_flat.tif")
summary(values(hs))
summary(values(hs)/cos((90.0-80.0)*pi/180.0))

I will try to parse the source code but it's been ten years since I've used C++.

Can someone with general knowledge about hillshades point me the right direction?

I thought I could get a "percent of the tile that is shaded" from running a hillshade. Upon further reading a hillshade outputs the illuminance of each tile. So, the cosine is only part of the equation, I would need to pull in the slope and aspect as well for additional trigonometric readjustments, and maybe more.

Turns out there are tools that are functionally equivalent for that I need called "shadow masks" that only output 1 or 0, indicating whether the tile is shaded.

This nail biting article referenced in gdal's source code helped enlighten me, in addition to this question on cloud shading, this question on solar maps, and this question on r.sun where a contributor indicates that r.sunmask is slow. A couple available tools:

GRASS r.sunmask

RSGISLib Elevation Module: shadowmask

USGS ISIS Shadow

Of course I need the speed of gdal's hillshading so I'll keep this question up until I or others come up with a solution that can process a 1000 square mile area in a second.


texts: Openlayers2 and Geoserver: OSM/Google Maps and WMS Overlay

texts:
I am experiencing projection problems with the following set up:

GeoServer with ESPG:4326 based shapefiles (*.shp)
OpenJump to edit my shapefiles
Openlayers with EPSG:900913 base layers
GeoServer provides the shapefiles as WMS to OpenLayers.

Here is my source code:

    <title>Openplanning mit Geoserver</title>
    <meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">

It seems that GeoServer is converting from ESPG:4326 to EPSG:900913 internally. Here you can see a sample request for Vienna (which is sent by my OpenLayers application to GeoServer):

http://localhost:8888/geoserver/wms?LAYERS=topp:ma_vienna&SRS=EPSG:900913&TRANSPARENT=TRUE&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&STYLES=&FORMAT=image/png&BBOX=1800244.8899219,6116147.6705496,1839380.6483984,6155283.4290261&WIDTH=256&HEIGHT=256

My problem: the requested file looks distored, as you can see here:

In case I open the WMS layer in the OpenLayers preview of GeoServer, it looks OK:

The question is: Why is my shapefile distorted in my OpenLayers application, but not in the GeoServer preview? As you can see, the GeoServer preview uses EPSG:4326 and my OpenLayers application EPSG:900913 (Update: EPSG:4326 as well; mouse coordinates are displayed as latitude/longitude). But, however, as GeoServer provides a WMS interface with EPSG:900913 capabilities, shouldn't it correctly reproject my vector sources?


texts: QGIS 3.14 vs "no Qt platform plugin could be initialized"

texts:
Ever since updating to the most recently released QGIS stable version, 3.14.15, I am getting the error pop-up window:
"This application failed to start because no Qt platform plugin could be initialized..."
I have doublechecked the contents of the Qt platform pathway (C:\OSGeo4W64\apps\Qt5\plugins\platforms) and do have all 4 DLLs. After hitting a brick wall I eventually performed a complete uninstall and reinstall of QGIS 3.14 but am still unable to open QGIS with Python for either stable or nightly releases.
Any help would be great.


texts: Need to fill the area of shapefile that contains LineStrings, which form a polygon, in ipython

texts:
The porblem I have is that I fetched the data about the administrative level of vienna (from OSM via overpass --> vienna relation id = 109166) as xml format. (json didn't work it had encoding problems) As I did not find any information of how to get it directly, I had to see if I could create the filled polygon with the data provided. With the ways included, forming a polygon itself doesn't work as some ways only have two nodes. Therefore, I created LineStrings in order to generate the border of vienna and saved it as a shapefile.

It looks like this:

Indicated by the different colors, one can see the LineStrings that form a polygon.
Now I would like to fill the area so that I can use it as a filled polygon. I am using ipython (with python 2.7) and currently searching for answers for the best approach to do this. Would be nice if there is a function that ipython can interpret all LineStrings and see that they are closed so it can write it as a polygon. Then I could identify it as polygon in order to fill the area. Working with programs like QGIS is not an option as I want to avoid manual corrections.

Dropbox link for shapefile

What I want is just a filled polygon of this relation in ipython.


texts: Software for analyzing GPS tracklogs to extract human-readable placenames

texts:
I have a pile of GPX-formatted tracklogs for many countries where I've travelled over the years (the GPS-device has been continuously on and recording for most of the trips). I would now like to analyse the logs automatically to extract a list human-readable names for the places that I've visited. For example:

detect and name border crossings
detect and name inhabited places where I've stayed longer

Does there exist ready-made software (preferably free and open source) that can be used to perform such analysis?

(By "tracklog" I mean: a sequence of timestamped longitude-latitude pairs.)


texts: Where does the default Earth radius in ST_Distance_Sphere come from?

texts:
MySQL says in the docs for ST_Distance_Sphere

  Calculations use a spherical earth and a configurable radius. The optional radius argument should be given in meters. If omitted, the default radius is 6,370,986 meters. If the radius argument is present but not positive, an ER_WRONG_ARGUMENTS error occurs.

PostGIS says in the docs of ST_Distance_Sphere, (though the docs aren't accurate anymore)

  Uses a spherical earth and radius of 6370986 meters.

Where did they get the default 6,370,986 meters from? WGS84 says major-axis radius is 6,378,137.0 m. PostGIS which now uses a Average Radius essentially uses 6371008.

Looking at the code

#define WGS84_MAJOR_AXIS 6378137.0
#define WGS84_INVERSE_FLATTENING 298.257223563
#define WGS84_MINOR_AXIS (WGS84_MAJOR_AXIS - WGS84_MAJOR_AXIS / WGS84_INVERSE_FLATTENING)
#define WGS84_RADIUS ((2.0 * WGS84_MAJOR_AXIS + WGS84_MINOR_AXIS ) / 3.0)

that means

-- SELECT 6378137.0 - 6378137.0 / 298.257223563;
WGS84_MINOR_AXIS = 6356752.314245179498
-- SELECT ( 2.0 * 6378137.0 + ( 6378137.0 - 6378137.0 / 298.257223563) ) / 3.0;
WGS84_RADIUS = 6371008.771415059833

Newer versions are much less efficient, more complex, and use Pro4j but they seem to do the same thing.

Still where does 6370986 come from?


texts: Existing standard country codes explicit in time and space

texts:
There are several comprehensive systems for coding countries such as ISO, TLD, M49. I am looking for a system that is explicit about the time period that the code refers to.
My application is a research database for time series modelling on a global scale. For example I need to put one data point in for Germany in 1989 and one for Germany in 1990. If i use the ISO standard codes DE refers to all of Germany after 1990, but only west Germany before 1990.
Does anyone know of a system that not only updates for political changes, but includes a record of the previous iterations?

Update: the UN statistical division M49 three number code has independent codes for political entities before and after major changes. I'm still looking for something with clear ties to geography.


texts: Google Earth Engine. Calculate Cloud and Cloud shadow for a especific region in am image collection of LANDSAT 8

texts:
I am working with hydrosedimentological data from a lake in my region. we are generating data on the concentration of suspended sediments, from the red band, for the entire history series of landsat 8 images. in this sense we would like to generate the percente of cloud cover data for the specific region of our lake, for a better selection of the data to be analyzed. I'm having trouble calculating the pixel_qa area of water for the lake and adding this data to the properties of the images in the collection,can someone help me?
this is my code:

 // region of Guaiba lake
 var region_guaiba = ee.FeatureCollection('users/eduardonmanara/Reg_guaiba');
  Map.centerObject(regiao_guaiba,10);

// data range
var start_date = ee.Date.fromYMD(2013, 1, 1);
var end_date = ee.Date.fromYMD(2020, 12, 31);

 // LANDSAT 8
var col_landsat_8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
  .filterDate(start_date, end_date)
  .select([1,2,3,4,5,6,7,10], ['BB','BG','BR','BNIR','BSWIR_1','BSWIR_2','BT','pixel_qa'])
  .filterBounds(regiao_guaiba);

// helper function to extract the QA bits
function getQABits(image, start, end, mascara) {
 // Compute the bits we need to extract.
 var pattern = 0;
 for (var i = start; i <= end; i++) {
 pattern += Math.pow(2, i);
 }
 // Return a single band image of the extracted QA bits, giving the     band a new name.
 return image.select([0], [mascara])
 .bitwiseAnd(pattern)
 .rightShift(start);
}

var col_landsat_8 = col_landsat_8.map(maskQuality);

// calculate the percentage of clouds inside the Guaíba lake

var getPercentCCG = function(image){
   var QA = image.select('pixel_qa');
   var water = getQABits(QA,2,2,  'water');
   var area_water = water.multiply(ee.Image.pixelArea()).divide(1e+04);
   var regiao_guaiba = ee.FeatureCollection('users/eduardonmanara/Reg_guaiba');
   var sum_water = ee.Number(area_water.reduceRegion({
  reducer:ee.Reducer.sum(),
  geometry: regiao_guaiba,
  scale: 30,
  maxPixels:1e13
   }))
   var guaiba_area = ee.Number(480.9)
   var percente_cloud_guaiba = sum_water.divide(guaiba_area);
   return ee.Feature(null, {'CLOUD_CLOVER_GUAIBA':percente_cloud_guaiba})
}

var col_areas = col_landsat_8.map(getPercentCCG);
print('coleção de áreas:',col_areas);
}

but the code dont work. I need to add another propertie for images in imageCollection containing the percentage of covering the lake itself.
then apply the cloud filter in an imageCollection contaning cloud cover las than 50% over the lake.
// A function to mask out cloudy pixels.
var maskQuality = function(image) {
 // Select the QA band.
 var QA = image.select('pixel_qa');
// Get the internal_cloud_algorithm_flag bit.
 var sombra = getQABits(QA,3,3,  'cloud_shadow');
 var nubes = getQABits(QA,5,5,  'cloud');
 //var cloud_confidence = getQABits(QA,6,7,  'cloud_confidence');
 var cirrus_detected = getQABits(QA,9,9,  'cirrus_detected');
 //var cirrus_detected2 = getQABits(QA,8,8,  'cirrus_detected2')
 // Return an image masking out cloudy areas.
 return image.updateMask(sombra.eq(0)).updateMask(nubes.eq(0).updateMask(cirrus_detected.eq(0))
 );
}


texts:
I search a lot about this problem but I don't find a solution, so...
I'm using a projection and the scaleLine display the value 100.000 greater than the real value, How can I solve this problem?

Here my OpenLayers 4 code:

Proj4js.defs["EPSG:27700"] = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs";

var projection = new ol.proj.Projection({
code: 'EPSG:27700',
units: 'm'
});

controls: ol.control.defaults({
  attributionOptions: ({
  collapsible: false
  })
}).extend([
  // add the actual scale in "m"
  new ol.control.ScaleLine({units:'metric', target:
document.getElementById('map_extras')})
]),
layers: [layers, orders],
target: 'mapdiv',
view: new ol.View({
projection: projection,
    center: ol.proj.fromLonLat([parseInt(setPos[0]), parseInt(setPos[1])], projection),
resolution: 2.8,
resolutions: [2.8, 1.4, 0.7, 0.28, 0.14, 0.07],
})

Also, something that I saw is when you move the map in x y axis the scaleLine sometimes changes without modifying the actual zoom, is this also a common problem when you are using a projection?


texts:
I am working on GeoServer 2.13 version. When I am working with WFS-T using OpenLayers I am getting WFS version 1.1.0 request, but my GeoServer accepts only WFS 2.0 requests.

How should I change the WFS default version to 1.1.0 in GeoServer?


texts: Trying to apply inlabru:::bru_fill_missing to a SpatialPointsDataFrame to fill in NAs

texts:
I'm struggling with a function from the inlabru R package, which could help me fill in the NAs in my data.
According to the helpfile, I need:
bru_fill_missing(
  data,
  where,
  values)

data is a SpatialPointsDataFrame, which I have - it contains a data.frame of 17 climate variables, some of these include NAs.
where should be "A matrix, data.frame, or SpatialPoints or SpatialPointsDataFrame, containing the locations of the evaluated values" - I'm guessing this would be the coordinates of the SpatialPointsDataFrame?
values  is "A vector of values to be filled in where is.na(values) is TRUE" - so this should be one of the climate variables from the SpatialPointsDataFrame.
cec2 <- inlabru:::bru_fill_missing(grid2, grid2@coords, grid2$cec)

...gives me:
Error in inlabru:::bru_fill_missing(grid2, grid2@coords, grid2$cec) :
  inherits(data, c("SpatialPixelsDataFrame", "SpatialGridDataFrame")) is not TRUE

Why is this not working? I seem to be using the right sort of objects:
> class(grid2@coords)
[1] "matrix" "array"
> class(grid2)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
> class(grid2$cec)
[1] "numeric"
> grid2
class       : SpatialPointsDataFrame
features    : 12179
extent      : 499974.9, 1928094, 8696615, 9781413  (xmin, xmax, ymin, ymax)
crs         : NA
variables   : 17
names       :     gbifID,             Species,    ID, Latitude, Longitude,    costpath,         cec,      Temp10,  thermicity,         tri,        fpar,        lai,        Elev,   soc_0.5cm,         sq3, ...
min values  :         49, Abelmoschus manihot,     1,  -11.521,       141, 0.158908026, 4.099550837,           0, 12.01840072, 0.160710489, 5.959515748, 1.23322692,           0, 43.46748724,           0, ...
max values  : 3039118788, Zygogynum argenteum, 13690,  -1.9667,   153.985, 422.5197024, 426.9292313, 101.4610783, 726.4499222, 3756.437675,         254,        254, 4204.876009,  1324.01969, 5.413501295, ...


texts:
I have run into a problem installing pgRouting. Following Anita Grasers guide, I have gotten as far as to execute the three .sql files.

The first sql file (routing_core) runs fine, but for the second and third (routing_core_wrappers and routing_topology), I get the following error message:

ERROR: type "geometry" does not exist
SQL state: 42704

I am using PostGreSQL 8.4 and have tested using PostGIS 1.4.2, 1.5, and 2.0. At first, I used a newer version of PostGreSQL, which ran into trouble at an earlier stage (thanks to all the answers provided to other posts, which aided me in getting me as far as this)

Any advice would be greatly appreciated!


texts: how do I get percentages of classes in a classified image?

texts:
I am using ArcMap 10.3

I have two layers, one classified raster, Landsat based dNBR, where classes show burn severity, and a shapefile containing polygons, drawn based on the dNBR image and that shows the extent of the burned areas (to have the area)

I need the percentages of the burn severity classes inside the polygons.
for example:
- low burn = X%
- medium = Y%
- high = Z%

but considering that the total area is the area of the polygons and not the whole image.


texts: Is there a way to add turn restrictions in A* and Dijkstra?

texts:
We are currently working with pgRouting and found that there's no way of implementing turn restrictions (e.g. banned right or left turns). While it is possible to assign "to_cost" and "rule" in Shooting * algorithm... I couldn't find a way of implementing these restrictions using A star and Dijktra algorithms.

Is there a way of implementing special turn restrictions in A star and Dijkstra methods ?


texts: Join between census shapefile and census CSV yields attribute table columns populated in disorderly fashion

texts:
I am using QGIS 3.10 to display American Community Survey variables at the census tract level. Last week, I successfully joined three ACS tables to duplicated copies of census tract geometries (same source), styled them, printed maps using the style, saved the project, and closed the application.
Today, I reopened the project to find that the original attributes from the joined attribute tables are no longer in the expected fields, and most geometries were initially not displayed in the map viewer. However, two attributes derived using field calculator remained in a clean row/column format, and the geometry seems ok, because all expected shapes were shown when reclassifying by one of the derived variables.
Both files are encoded in UTF-8, and no changes were made to the files between the last time I used them and present. Several other duplicate copies of the shapefile joined to different census tables were not affected. Is it apparent from this description what happened here, and if so, how to keep it from happening in the future?


texts: QGIS Join Attributes by Location not working from Python script

texts:
I have a shapefile with points, and a shapefile with polygons. I want to sum the attribute value of all points that are within a certain polygon, and I'd like to do this for each polygon in my dataset. QGIS Join Attributes by Location tool is helping me to achieve this, and is working perfectly when using the tool via the Processing Toolbox. However, I would like to include the tool in an automated script. I tried many different arguments, but my code is (still) not working, and results in a dataset where the new polygon attribute has NULL values.

I'm using the following code, which should be exactly the same as the tool is running when I'm accessing it from the Processing Toolbox (I even tried copying and running the exact same code from the History & Log window):

sum = processing.runalg("qgis:joinattributesbylocation","polygonlayer.shp",
                        "pointlayer.shp","['contains']",1,"sum",1,None)

I'm using QGIS 2.8.2 Wien.


texts: Adding time as 'year' band to image collection

texts:
I have uploaded 35 images (annual maximum NDVI) to GEE and put them manually (drag and drop) together to an image collection. There is no time band, but for the ee.Reducer.linearFit() function, I need to add the year of each image as a time band to each image.
I found this image.addBands(image.metadata('system:time_start') but since there is no system:time_start, this can not be used.
This is another way of doing it:
var Maxyears = ee.List.sequence(1981, 2015);
var NDVIMax35 = ee.ImageCollection.fromImages(
  ee.List.sequence(0, NDVI_Max_Stack.bandNames().length().subtract(1))
    .map(function(i) {
      return NDVI_Max_Stack.select(ee.Number(i))
        .rename('Max35')
        .set('year', Maxyears.get(i));
    })
  );

But I have the image collection already, so I need a different function than fromImages()
https://code.earthengine.google.com/?scriptPath=users%2Fwelp1992%2FMaster%3AClimate%20Data


texts: Setting visibility of columns of Attribute Table in QGIS?

texts:
Is there a possibility in QGIS setting a field (column) of the Attribute Table to not visible?

In the properties-dialogue of the layer, I couldn't find anything.


texts: Select features from attribute table and then create a new field with the selected attribute data

texts: I have a shape file which contains land use typology for my administrative area. I’m investigating deficiency in certain land uses and there are 2 ‘standards’ that we have to use. These standards use different land uses. The Shapefile has a field in the attribute table that contains all the land uses, for example each polygon can have a land use of either A1 – 12, B1 – 12, C1-12 etc. One standard  requires land uses A1 – 6 and C5 -9 and the other A1-12 and B1-9. The Land uses A, B & C have separate fields. I have created an SQL query in the advanced search that selects these land uses. What I now want to do is create 2 new fields, one which populates the selected land use for each standard. I know I can select the attributes and save as a new selection but ideally I would like to populate 2 new fields with the attributes from the existing land use field.


texts: Which area measurement to trust? - The same polygon in GoogleEarth v. ArcGIS

texts:
Possible Duplicate:
  What is the most accurate coordinate system for calculating areas of polygons?

I'm currently involved in comparing the "official metropolitan" areas of several cities against their actual overall urban footprint. The official metro areas' areas can be found in other sources, but I need to get the most recent urban footprint area.

I've found the fastest way to do this is to draw a polygon in Google Earth around a city's urban sprawl. As I want to create an explanatory graphic, I then brought my GE polygon into ArcGIS 10, where I converted it from KMZ to layer, and then saved the layer data as a new feature class in my specific city/map geodatabase.

Now, out of habit, b/c I'm focusing on a specific small area, I'm projecting things in UTM. I checked the attribute of my new "urban footprint" feature class polygon and found there to be an ~600,000m2 difference between the polygon area in Arc and the same polygon's area given in Google Earth.

In light of this, my question is which area value is more accurate and should I use in future calculations I'm going to do?


texts: QGIS crashes when setting new Symbol to layer with python

texts:
I followed the QGIS 3 cookbook and created a new custom symbol layer class. But when I try to create a new the symbol, QGIS crashes. I create it like this.

symbol = QgsMarkerSymbol({FooMarker(20)})

The link for the cookbook is the following: https://docs.qgis.org/testing/pdf/en/QGIS-testing-PyQGISDeveloperCookbook-en.pdf page: 31

My FooMarker class looks like this.

class FooMarker(QgsMarkerSymbolLayer):
    def __init__(self, radius=4.0):
        QgsMarkerSymbolLayer.__init__(self)
        self.radius = radius
        self.color = QColor(255,0,0)

    def layerType(self):
        return "FooMarker"

    def properties(self):
        return {"radius": str(self.radius)}

    def startRender(self, context):
        pass

    def stopRender(self, context):
        pass

   def renderPoint(self, point, context):
        color = context.selectionColor() if context.selected() else self.color
        p = context.renderContext().painter()
        p.setPen(color)
        p.drawEllipse(point, self.radius, self.radius)

    def clone(self):
        return FooMarker(self.radius)


texts: Measuring true length of polyline over terrain in QGIS?

texts:
I use QGIS for a variety of mapping purposes and also for managing UAV/drone data. I need to calculate the true length of a mapped fence over varying terrain for which i have a detailed DEM.

Is this possible?


texts: Wrap CSV lat/long points around 180 meridian in QGIS

texts:
I have a CSV file with lat/long coordinates. I also have a custom CRS to center the map on the Indian Ocean. Is there a way for me to force the points to plot on the blank band of the map where the 180 meridian is?

My proj4 string is +proj=eqc +lon_0=90 +datum=WGS84 +units=m +no_defs +lon_wrap=90. +lon_wrap=90 doesn't appear to have done anything.

As a temporary solution, I first converted my CSV points to a points shapefile, then rasterized my shapefile at the appropriate resolution. In my case my CSV file contained lat/long coordinates every 1/2 degree, so I set my resolution to be 0.5 map units = 1 pixel, because my projection was in degrees lat/long. My raster wraps with no issues:

However, as another user noted, this is only a visual solution and doesn't actually solve my original problem.


texts: Counting number of features within another feature using aggregating function in field calculator of QGIS?

texts:
My data: 1 polygon layer + 1 point layer
My target: Calculate the number of points and the sum of a field within each polygon using the aggregate function in the QGIS field calculator.

I know I can simply use the QGIS functions Join attributes by location (summary) or Count points in polygon but I´m interested in the aggregate function in the field calculator since I don´t understand exactly its meaning.

Using

aggregate('hh', 'sum', "sum_pop")

'hh' - the name of the point layer
'sum_pop' - the name of the attribute field from the hh layer (sum of the population of a point)

I get the sum of the whole field in the attribute table, so each feature has the same value.

How can I change the code for the aggregate function to calculate the sum of a value based on the number of points in each polygon feature and the number of features within a polygon so I get different values for each polygon?


texts: Combining two aligned overlapping rasters in raster calculator

texts:
I have two rasters with the same extent, projected coordinate system, cell size and they are aligned.

"raster1" has continuous values ranging from 0 to 6.09
"raster2" has continuous values from 0 to 91.25

What I need to do is creating a new raster in which, cells that have "value=0", based on "raster1", they will remain 0 but in any other case, cells will be assigned a value based on the sum of the values from "raster1" & "raster2".

The expression I use in Raster Calculator is :

Con(("raster1" == 0), 0, (Con(("raster1" > 0) & ("raster2"  >= 0)), ...?

Obviously, this is not complete.  What should I change?


texts: Generate Raster Mosaics in a folder using all TIFFs in the folder, then iterate in the next folder

texts:
I have 200+ directories, each directory contains x amount of TIFFs which I want to mosaic.
I have a snip of code that will list the TIFF contents of all dirs under a parent directory, now I just need to parse the list (per dir), pass it to the mosaic program and run all TIFFs in the first dir, then loop to the next dir in the list and repeat the mosaic operation.
import arcpy
import os

#Set the workspace enviorment setting
arcpy.env.workspace = r"C:/Work/SpaceNet/Test2"

walk = arcpy.da.Walk(datatype="RasterDataset")

for dir_path, dir_names, file_names in walk:
    for filename in file_names:
        print(os.path.join(dir_path, filename))

C:\Work\SpaceNet\Test2\00114 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190822100114_20190822100413_tile_678.tif
C:\Work\SpaceNet\Test2\00114 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190822100114_20190822100413_tile_679.tif
C:\Work\SpaceNet\Test2\00114 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190822100114_20190822100413_tile_680.tif
C:\Work\SpaceNet\Test2\00628 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190823100628_20190823100943_tile_5430.tif
C:\Work\SpaceNet\Test2\00628 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190823100628_20190823100943_tile_5431.tif
C:\Work\SpaceNet\Test2\00628 SN6_Train_AOI_11_Rotterdam_SAR-Intensity_20190823100628_20190823100943_tile_5432.tif


texts: How to handle multiple layers when importing a polygon shapefile into GRASS?

texts:
When I import a shapefile into GRASS, it seems to arrive as 3 layers.
I.e. they appear in the qGIS toolbox browser as:

myMap 0_polygon, myMap 1_polygon, myMap 2_polygon

Is there a way to merge the layers into one?
I have run:

v.clean tool=snap,break,rmdupl,bpol thresh=10 input=myMap output=myMapCleaned

And searching on this topic suggests that I should now run v.patch to merge the layers. But how do I specify the individual layers in the 'input=...' string?
Thanks.


texts: Different maps/views/lay-out in 1 project

texts:
I am looking for a option in QGIS to make different maps/views/layout in one project. I have made a map in AutoCAD of houses, roads and public space. I've set all the object in different layers so I can turn the objects on and off in QGIS.

I want to make different views with specific layers turned on/off without interacting with other views.

Lets say 1 map with only the houses, 1 with only the roads and 1 with only the public space in the same project.


texts:
I am attempting to contour a DEM raster using gdal.ContourGenerate().

This is my code:

in_data_set = gdal.Open(TEMP_CROP_DEM_PATH, gdal.GA_ReadOnly)
in1 = in_data_set.GetRasterBand(1)

ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(TEMP_CONTOURS_PATH)
contour_shp = ogr_ds.CreateLayer('contour', dest_srs)

field_defn = ogr.FieldDefn("ID", ogr.OFTInteger)
contour_shp.CreateField(field_defn)
field_defn = ogr.FieldDefn("elev", ogr.OFTReal)
contour_shp.CreateField(field_defn)

gdal.ContourGenerate(in1, contour_offset, 0, [], 0, 0, contour_shp, 0, 1)

This code has worked for me on several tests, however it does not work when there is an area of no data in the DEM. When contouring a DEM with no data areas, QGIS freezes on the last line (ContourGenerate) and crashes.

I have looked through the GDAL API and made sure all my arguments are correct. I have tried changing the 5th argument useNoData to 1 and saw no change.

Running in1.GetNoDataValue() also causes QGIS to crash.

I am running out of ideas of where to look to solve my problem. Has anyone come across something similar? Are there other ways to achieve my goal, create vector contours from a DEM in Python?


texts: Calculating the shortest distance coordinate to vector layer

texts:
I want to calculate the nearest distance from each coordinate in a vector layer containing more than 17 000 coordinates to another vector layer. I tried the GRASS v.distance tool but it gives the following error:

Is there any other method to do this?


texts: Generating choice list (dropdown) in toolbox to allow user to select field in CSV as input?

texts:
I need to input a csv as parameter 1. Parameter 2 should be a dropdown displaying the various fields (the headerline) where I can select a column to use for geoprocessing.

Param1:  c:\my.csv
Param2: (dropdown) select header


texts: Is it right to remove salt & pepper effect from a product produced by clean images?

texts: I have two images from classification (different years) and I want to find the changes of land cover . I have cleaned the two images (salt and pepper). When I create  a cross product of the category values from these images I get an image with a lot of "noise". I want to know If it is right to clean that image or not(because it is a product of two clean images). I use grass and the final image created with "cross".


texts: Extract points from regular data.frame with coordinates based on SpatialPolygonsDataFrame in R

texts:
I have a regular data.frame called datos.merged

  class(datosmerged)

[1] "data.frame"

The last two columns are the corresponding geographical coordinates

  datosmerged[1,c(23,24)]

   XCOORD   YCOORD
1 2473520   668648.7

I have a file called shape

  class(shape)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

I would like to use this very simple shape to subset the previous data.frame. How can this be done?


texts:
I'm trying to write a query that excludes all entries with the zone AG but would include any blank entries. I've tried the following but it doesn't work:

Zone <> "AG" and Zone = "Null"

Any suggestions for how to build this?


texts:
While I was travelling with my GPS "ON" and tracking my route I found this strange lines in pink and green dotted lines when I zoomed out the map I found that the lines where coinciding with the airport(RGIA AIRPORT HYDERABAD,INDIA) runway lines which is just around 5km from my location.

So, my question is if, this are the runway approach lines then why would it indicate or show in google maps.

Co-ordinates: 17°14'10.5"N 78°28'23.9"E


texts: Why is the result of a raster calculation integer when it should be decimal?

texts:
I'm trying to calculate the NDVI using raster calculator, but the result is integer instead of decimal.

NDVI = (NIR - RED)/(NIR + RED)

The result is normally in the range -1 to 1, but when i tried to calculate de result is -1, 0 and 1, and it does not have have any number in between like, -0.9 or 0.1. And i need those values so i can classify the result in term of soil occupation


texts: Converting a KMZ file containing an api link for use in python's geopandas

texts:
I've got a KMZ file which contains an API link for a map overlay, this is the link to that file. (The corresponding json file (link here), posted  only gives me a link to a low resolution image and the bounds.)

I want to use this map in some python application, mainly using geopandas, which uses Fiona which does not support KMZ.

I've found some web applications that can convert KMZ files to geoJSON or geoTIFF, but not this KMZ file (or I get 1 pixel or something).

So I tried installing Google Earth Pro and loading the file into that, which works just fine, it seems Google Earth knows how to interface with this. However, when I then try to export the layer I end up with the exact same KMZ file I started with.

Is there anyone who can point me towards a program or script that can do this conversion for me?


texts: Connecting to Excel spreadsheet through ArcPy giving ERROR 000582

texts:
My ultimate goal is to automate kernel density analysis on xy data from a table. I am trying to learn Python but I am very new. I built a simple model to take table data and display xy data. When I attempt to export the model as a Python script, the script fails. I do not have a programming background but I have programmers available to me.

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# test.py
# Created on: 2014-02-12 13:14:37.00000
#   (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy

# Local variables:
Input_Table = "C:\\Users\\cbriglia\\Desktop\\GIS\\Worksheets\\accidents_jan_2013.xls\\Sheet1$"
Sheet1__Layer = "Sheet1$_Layer"
Ouput_Feature_Class = "C:\\Users\\cbriglia\\Desktop\\GIS\\Shapefiles\\acc_jan_2013_builder.shp"

# Process: Make XY Event Layer
arcpy.MakeXYEventLayer_management(Input_Table, "longitude", "latitude", Sheet1__Layer, "GEOGCS['GCS_North_American_1983_CSRS',DATUM['D_North_American_1983_CSRS',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision", "")

# Process: Copy Features
arcpy.CopyFeatures_management(Sheet1__Layer, Ouput_Feature_Class, "", "0", "0", "0")

Runtime error  Traceback (most recent call last):   File "", line 19, in    File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 6348, in MakeXYEventLayer     raise e ExecuteError: ERROR 000582: Error occurred during execution.


texts: Customizing map.on("click") while using default editor widget

texts:
My question is whether or not it is possible to use the default editor widget in my case given the needs of my application. My map consists of points displayed from four different FeatureLayers. They represent four different types of traffic counts, and my application allows users to associate/dissociate them with traffic count data requests. (Count data requests are non-spatial records stored in a data table.) They are displayed on my page below the map in a Dojo DataGrid.

If no traffic count request is selected in the DataGrid, the Editor and TemplatePicker should behave as normal (i.e. they should select points, display the Attribute Inspector, and allow the points to be dragged around). If a request is selected, the Editor should be temporarily disabled so that clicking the point toggles it as being associated/dissociated with the traffic count request clicked in the DataGrid. The AttributeInspector should not be displayed

As it stands right now, if there is a request selected in the DataGrid, clicking a point on the map still brings up the AttributeInspector. Is it possible to customize the Editor in this way?

The definition of my TemplatePicker and the logic of what I want map.on("click") to do is below.

Thanks

function generateTemplatePicker(layers) {
    widget = new TemplatePicker({
        featureLayers: layers,
        rows: "auto",
        columns: 2,
        grouping: false,
        style: "width:98%;"
    }, "templatePickerDiv");

    widget.startup();

    widget.on("selection-change", function() {
        unselectRequest();
        selected = widget.getSelected();
    });
}

map.on("click", function(evt) {
        if (selected) {
            var newPoint = new Graphic(evt.mapPoint, selected.symbol, null);
            selected.featureLayer.applyEdits([newPoint], null, null);
            widget.clearSelection();
        } else {
            if (!evt.graphic) {
                if (requestEditingMode) {
                    unselectRequest();
                }
                else {
                    //*HIDE ATTRIBUTE INSPECTOR*//
                }
            } else if (!requestEditingMode) {
                //*SHOW ATTRIBUTE INSPECTOR*//
            } else {
                var pointAlreadyRequested = pointIsRequested(evt.graphic._graphicsLayer.id, evt.graphic.attributes.OBJECTID);
                if (pointAlreadyRequested) { removePointFromRequest(evt); }
                else { addPointToRequest(evt); }
            }
        }
    });


texts:
I would like to take a ArcMap mapbook and export each individual page to a PDF named for the page and/or GridID.  I've never messed with Python...  But I see that the arcpy.mapping module should do the trick.  Doesn't seem like this should be too difficult, but I am failing...

i took Gistech007's code and modified by adding a line  -
pageName = mxd.dataDrivenPages.pageRow.stand_key
and then instead of having str(pageNum) in the name I have str(pageName).  so it looks like this now:

import arcpy
 print "Setting Map Document. . ."
  mxd = arcpy.mapping.MapDocument("CURRENT")
  for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
  mxd.dataDrivenPages.currentPageID = pageNum
  pageName = mxd.dataDrivenPages.pageRow.index
  print "Exporting page {0} of {1}".format(str(mxd.dataDrivenPages.currentPageID), str (mxd.dataDrivenPages.pageCount))
  arcpy.mapping.ExportToPDF(mxd, r"D:\Temp\Page_" + str(pageName) + ".pdf")
 del mxd

Where "index" in the line - pageName = mxd.dataDrivenPages.pageRow.index is a field in my table.  This runs and works just fine, but inst what i am after.  I need the PDF to be named with only the str(pageName) and no prefix of "Page_".  If I remove the "Page_" from the code it throws an error... How can i create the PDF and have it named using only the value from  pageName?


texts: File geodatabase not fully loading in QGIS?

texts:
I am very new to GIS. I have been trying to load a gdb file into QGIS, but one of the layers appears completely different from the description in the metadata. The data can be downloaded here. https://pubs.er.usgs.gov/publication/ofr20181088

The layer that doesn't work for me is called "Recent_2013_2017_Mining_Activity_Dataset". It's supposed to contain loads of points where mining activities take place in the Central African Republic, but the layer I get in QGIS is a grid consisting of tiny squares covering the whole country.

Would anyone mind having a look at this for me?

I really don't know what is causing this.


texts:
I want to change attribute table columnd value when intersection is found but this is not updating. see comments for details

mapcanvas = iface.mapCanvas()
layers = mapcanvas.layers()
print layers[0].name()
print layers[1].name()
selections = []
for f in layers[0].getFeatures():
  for a in layers[1].getFeatures():
    if a.geometry().intersects(f.geometry()):
      intersection = a.geometry().intersection(f.geometry())
      print 'polyTwo  '+a['PolyTwo']
      a['PolyTwo'] = 'xyz'#updating value PolyTwo to xyz
      print 'After polyTwo  '+a['PolyTwo'] #change value is showing here but not in db
      print intersection.exportToWkt()
      selections.append(f.id())


texts:
I've discovered, that I have no idea, how height is usually represented in GIS.

For example, let's say we have software, that lets you create a mission plan for an UAV such as

take off
go to point x1
go to point x2
...
land

Where "point" is represented by

latitude
longitude
height

Latitude and longitude are "intuitive" (there's no question about what they actually mean, once you select a geodetic datum). However, height is not. How can one represent it? Things, that come to mind:

height in relation to takeoff point
height in relation to "sea level"
height in relation to some other reference point
distance from the center of the Earth (I mean, why not?)

Is there some canonical way of representing heights, preferrably, a way that is not tied to the problem domain (UAVs, etc)?

Of course, I am ready to accept an answer, that will point me to some reference / literature on the topic, instead of answering the question directly. I am, in fact, quite ignorant about GIS, so I don't have much of an idea even which keywords to search for.


texts: Geopandas performance appears quite slow

texts:
I'd like to confirm the performance I am observing when utilizing distance from a Geopandas Series.

Problem

Time to complete operations appears to be far greater than comparable operation in PostGIS. I would like to understand if this is known and, if so, if there are suggestions as to methods for making Geopandas more performant, particular regarding geo operations like buffering.

Goal

Given a set of geometries, I would like to calculate some aggregate (e.g. sum()) for each geometry plus those that fall within a given distance of the reference geometry.

Example

Given a set of geometries, I would like to calculate a quarter mile buffer (402 meters) around them and gather the sum the attribute du (dwelling units).

Current Strategy

Current method, utilizing solely the centroids in an attempt to be performant:

# precompute the centroids
centroids = df['wkb_geometry'].centroid

def test_measure(row):
    center = row['wkb_geometry'].centroid
    return df.loc[centroids.distance(center) < 402, 'du'].sum()

df.apply(lambda row: test_measure(row), axis=1)

Precomputing centroid and then using distance does appear to introduce some efficiencies compared to buffering and using within operation. Cost, regardless of methodology, grows at rate n^2 due to the fact that test_measure runs n times where n is the row count and must be run once for each of the n rows.

Times:

100 rows: 0.17s
1000 rows: 11.11s

Prior strategy

Prior, with buffering and within it took about:

100 rows: 1s
1000 rows: 100s

Prior method, not using centroids:

def test_measure(row):
    buffer_shape = row['wkb_geometry'].buffer(500)
    return df.loc[df['wkb_geometry'].within(buffer_shape), 'du'].sum()

Thoughts

Were I to perform a similar operation (using st_dwithin) in PostGIS, I would be able to run the operation in the following times.

Times:

100 rows: not run
1000 rows: not run
24000 rows: 75s

For reference, here is an example of that sort of SQL query:

CREATE OR REPLACE FUNCTION agg_within_dist(
    in_id int,
    in_geometry geometry,
    OUT id int,
    OUT du float)
AS
$$
    SELECT
        $1 AS geography_id,
        SUM(CAST(ref.du AS float)) AS du
    FROM s1.s1_scenario_final AS REF WHERE st_dwithin($2, ref.wkb_geometry, 402);
$$
COST 10000
LANGUAGE SQL STABLE strict;

SELECT (f).*
FROM (
        SELECT agg_within_dist(geography_id, wkb_geometry) AS f
        FROM scenario
     ) s


texts: How to find inner points near to corner of rectangular GeoJSON feature?

texts:
This question is in continuation of How to check if GeoJSON feature is in rectangular shape and find corner points?

I have a collection of GeoJSON features as Polygon & MultiPolygon which are saved in MongoDB. Many among them are in rectangular or square shape, while others are odd shaped. From above question's answer, I am able to find if the feature is rectangular or not also can find the corner points of it. But I actually need inner points which are near to corner (see attached image). I want to place label at that geolocation using Google map label library. How can I do that?

I tried to +/- some value from corner points but it is not working because some shapes are tilted towards left and some towards right.


texts: How to change the default color used by openlayers hover?

texts:
It is possible to change the color of highlightCtrl when I use onselectedHover? Default is light blue.

I want to change it to another one eg. red? Does anyone know where this value is defined?


texts: OpenLayers displays nothing when I try to add WFS Layer from Geoserver - can't work out where the problem is occurring!

texts:
I hope someone can help me with this problem as I'm very confused!

I recently downloaded the IHO Sea Areas shape file from http://www.marineregions.org/downloads.php

One of my colleagues has uploaded the shape file to PostGIS and served it up using Geoserver. Because I want to display it over a google maps base layer in OpenLayers, we tried to set Geoserver to "Reproject Native to Declared" from EPSG:4326 (detected by Geoserver) to EPSG:3857 (Google projection). However, when trying to access the layer through Geoserver, I got a 404 The requested resource () is not available.

The only way that Geoserver would serve the file was if we set it to "Force Declared", at which point the WFS would return the data as requested.

So, I don't really understand why that's the case, but I got a WFS that seemed to work.

Next step (note: I have set up OpenLayers with a declared projection of EPSG:3857 on the map), I added the layer to OpenLayers using:

IHO = new OpenLayers.Layer.Vector("IHO Areas",{
            style: {
                    strokeColor: "green",
                    strokeWidth: 1,
                    fillColor: "green",
                    fillOpacity: "0.2",
                    cursor: "pointer"
                },
                projection: "EPSG:3857",
                strategies: [new OpenLayers.Strategy.Fixed()],
                protocol: new OpenLayers.Protocol.HTTP({
                    url: "http://*MyWebProxy*/IHO_Areas.php",
                    format: new OpenLayers.Format.GeoJSON()
                })
            });

    IHO.events.on({"loadend":function(){console.log(this);}});
    map.addLayers([IHO]);

Now, the "loadend" event successfully fires, the logged IHO layer shows that it, and all its features, have been fetched successfully and the projection is logged as EPSG:3857 on the layer object, no errors are thrown by javascript/OpenLayers. But, nothing appears on the map.

This is especially weird as I copied and pasted the code from one of my earlier projects that worked perfectly using a different WFS vector layer source. So it seems there must be something weird about the IHO shape file or the data extracted from it, but I don't really know what to look for to fix it. Can anyone suggest anything to help?


texts: Pulling ArcSDE feature classes from current map to fieldmapping as DETable?

texts:
I'm working with SDE layers in arcmap. trying to avoid pointing to the sde file itself. I've tried workspace.Pathname, it was empty for sde remote workspace. I've also tried IName from workspace, empty again. Any ideas on what I should be doing differently? BTW, it works just fine for GDBs.

        Dim gputilities As ESRI.ArcGIS.Geoprocessing.IGPUtilities = New ESRI.ArcGIS.Geoprocessing.GPUtilitiesClass()
        'Dim inputTableA As IDETable = CType(gputilities.MakeDataElement(pWorkspace.PathName + "\Weld", Nothing, Nothing), IDETable)
        Dim inputtables As IArray = New ArrayClass()
        Dim gpValueTableObj As ESRI.ArcGIS.Geoprocessing.GpValueTableObject = New ESRI.ArcGIS.Geoprocessing.GpValueTableObject
        gpValueTableObj.SetColumns(1)
        Dim pdet As IDETable = New DETable
        For Each pLayer In checkedLayers
            Try
                gpValueTableObj.AddRow(pLayer)
                Dim workspaceFactory As SdeWorkspaceFactory = New SdeWorkspaceFactoryClass()
                Dim workspaceName As IWorkspaceName = workspaceFactory.Create(Nothing, pLayer, pWorkspace.ConnectionProperties, 0)
                Dim pName As IName = TryCast(workspaceName, IName)
                inputtables.Add(CType(gputilities.MakeDataElementFromNameObject(pName), IDETable))

            Catch ex As Exception
            End Try
            'inputtables.Add(CType(gputilities.MakeDataElement(weldLayer.Name, Nothing, Nothing), IDETable))
        Next
        fieldmapping = New ESRI.ArcGIS.Geoprocessing.GPFieldMappingClass()
        fieldmapping.Initialize(inputtables, Nothing)


texts: Show Geometric Network Topology data in the web

texts: I am trying to locate any API's for JavaScript (or other) that will help me to display current linear Geometric Network Topology and add multiple point/line geometry from JSON in HTML/JS. I specifically need to have my users perform typical tools like trace route and analysis in web interface on this linear network. I am just trying to find out if there is a JS or other Web API to handle these types of requests?


texts: Reassigning raster cell values from another attribute in raster attribute table?

texts:
I have a large (size of lower 48 U.S.) 100m resolution raster dataset that has a raster attribute table (RAT). The data were created by someone else. One of the attributes in the RAT has the data values I care about, while the actual raster cell values are meaningless. I've been trying to figure out how to reassign the value of the attribute I care about to the actual raster cell value but have not been able to figure out how to do this. The reason I want to do this is that I have a fairly involved custom Python tool that someone else developed that will do the calculations I need, but unfortunately it doesn't have the ability to use RAT attributes. It only uses the raster cell assigned value which as I said in this instance is meaningless.

I tried a roundabout way of running the raster to point tool to create a point feature class with the target field being the attribute I care about, thinking I could then create a new raster from the resulting point feature class. However, the raster to point tool failed stating I didn't have enough temporary memory (I'm running an Alienware machine with 16 GB of RAM and a 1 TB data drive).

Do I have any options? I'm using ArcGIS 10, but if I could do what I describe readily in another package I'd be willing to try another way (I'm not married to ArcGIS).


texts: Created features disappear when saving edits in ArcMap

texts:
I am trying to create point features in a template geodatabase. I am adding the features to a feature class in an edit session. When I save the edits all of the features I just created disappear and there are no attributes in the attribute table for that feature class.
I checked the properties, and the feature class I am editing has no extent but it is in a geographic coordinate system and the data frame is projected.
Could this be the problem and if so how do i fix it?


texts: Can I consider as continuous data vegetation layers on Maxent?

texts: I would like to know if we can consider as continuous data vegetal layers when disposed on seperate different .asc files (e.g.: Barren.asc; Water.asc; Evergreen Broadleaf Trees.asc; etc.)? The data on these maps seem to be in percentages (from 0 to 100) and I am working on Maxent. Depending on weither you consider the layers categorical or continuous, the output results are quite different. Also, it is mostly adviced to consider vegetation layers as categorical but I do not know if it is appropriate as the layers are in percentages and not categories. Here is a link to the layers I use: http://www.earthenv.org/landcover


texts: Why can't I add numeric(n,m) field values to PostGIS geometries in QGIS 1.7.1?

texts:
Actually I'm using QGIS 1.7.1 to create and visualise some geodata hosted by PostgreSQL 8.4 and PostGIS 1.4.0. The geodatatable is defined like this:

CREATE TABLE u2_baum
(
    baum_id integer NOT NULL,
    kronendurchmesser numeric(3,1),
    hoehe numeric(3,1),
    art_id integer,
    the_geom geometry,
    CONSTRAINT baum_id PRIMARY KEY (baum_id )
)
WITH (
    OIDS=FALSE
);

I added some datasets on database level, using some simple INSERTs.

Everything works fine, except adding attributes to new geometries in QGIS.  I am able to define the geometry and all integer attributes using the QGIS tools. The only problems were caused by the attributes with the numerical[n,m] fieldtype. One can't enter new values in these fields after creating a new geometry. The field can be "write activated" in the attribute table, but it isn't possible to enter values. Also it isn't possible to Copy and Paste some values (like 10.3) into this fields. But it is possible to change existing numerical[n,m] values.

Do you know this error and a way to fix it?

PS: Due to some restrictions I can't actually upgrade to QGIS 1.8


texts:
I have a series of plot points in the UK. All I want to do is overlay a simplified map of the UK, ideally just an outline.

I obtained some simplified colour maps from the OS website, added as a layer, but I can’t see the points beneath. I can’t merge them.  Perhaps the map has to be a vector rather than what I have , a raster.

I tried altering the layer order but no change. I am afraid.

If I "Zoom to layer" then it will show the layer of choice. But I can't see both together even if I alter opacity.

I think I am loading the map raster wrong, I am just loading it from the Browser and selecting the file directly. Perhaps it needs a bit more than that.


texts: Create VRT from multiple hdf files with Rasterio/GDAL

texts:
I have several Landsat 8 Scenes, each consisting of 11 bands. I want to mosaic them by creating a VRT and afterwards exporting the mosaic as TIFF again. I tried using GDAL via Python. But maybe it would be easier using Rasterio.

I tried:

rasterio.vrt.WarpedVRT(src_files_to_mosaic, "/.data/out.vrt", crs = "EPSG:32734")

But I get the error (one of several):

The WKT could not be parsed. OGR Error code 5

I think this should be an easy task, but it takes way to long me to figure out what to do.


texts:
I am new to GIS.
I'm trying to get a GeoJSON file for US ZIP3 borders to use with plotly.
I found a shapefile here and managed to simplify it using https://mapshaper.org/ and after repairing converted it to geojson using QGIS after repairing.
When using one or two of the ZIP3 polygons it seems to be working fine, but when using my full dataset  the whole map appears 'flooded'. Maybe one of the polygons is not terminated correctly?
How would one go about checking for this?


texts: Erode islands of pixels while polygonizing a raster

texts:
I am polygonizing a single band raster.

I am using gdal.Polygonize() function to vectorize the raster.

img= r"C:\Users\aoi2.tif"
src_ds = gdal.Open(img)
srcband = src_ds.GetRasterBand(1)

# create output datasource
dst_layername = "output"

drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername +".shp")
src = osr.SpatialReference(wkt=prj)

dst_layer = dst_ds.CreateLayer(dst_layername, srs = src)

gdal.Polygonize( srcband, srcband, dst_layer, -1, [], callback=None )

del src_ds, srcband, dst_ds, dst_layer

output shapefile looks like:

My aim is to erode the islands of pixels smaller than 3 pixels and keep the rest, how can I do that?


texts: Finding intersection topology errors in QGIS

texts:
I would like to find the intersection topology errors in a line layer using QGIS 2.18. In the topology checker I found no intersection errors. I know GRASS v.clean and v.build commands can solve these errors. I found that the Merge selected features tool from advanced editing toolbar (as a side effect) breaks intersecting lines and then using Multipart to singleparts from vector menu gives back the original single line features, but this way the attributes are lost.

After selecting all lines and merge them the intersecting lines are broken.

I would like to find a the intersection errors without automatically breaking the lines using QGIS. I need something similar to Topology checker.


texts: GeoTools: Calculate length along line from start vertex up to some point on the line?

texts:
I have single LineString and some Point (outside the line), I snapped it to the Line (using DistanceOp so as far as I understood the snapped point should lay on the line). How can I calculate length along the line from start vertex of the LineString up to the snapped point on the line?

    import com.vividsolutions.jts.geom.Coordinate;
    import com.vividsolutions.jts.geom.GeometryFactory;
    import com.vividsolutions.jts.geom.LineString;
    import com.vividsolutions.jts.geom.Point;
    import com.vividsolutions.jts.operation.distance.DistanceOp;

    GeometryFactory gf = JTSFactoryFinder.getGeometryFactory();
    Point outsidePoint = gf.createPoint(new Coordinate(outside.lat, outside.lon));
    Coordinate snappedPoint = DistanceOp.nearestPoints(outsidePoint, LineString)[1];

I though about this: LineString is a List of coordinates, Find the index of the vertex, closest to the snapped point, and replace it with the snapped point. Then construct a LineString from start vertex up to the index found. Then sum length calculated between each consecutive vertices of the line.

First problem is situation like this:

Second problem: When some sophisticated algorithm calculates length - Does inside happening something more complex than just sum lengths between 2 consecutive points?


texts:
I have a question
from a DXF file " topographic precision 3 cm " with 3d polylines to get one shp line with information 3d is possible?


texts:
I have an idea of what I want to do but I don't know yet what tools will bring me there.

I have administrative boundaries in one layer and coastline in another layer. The coastline is categorized in different classes. What I want to find out is how much of each category falls into each administrative boundary. So the outcome would be percentages of the administrative unit's coastline for each class. So say the coastline from unit 1 is 40% class A, 50% class B and 10%  class C. Unfortunately, the two shapefiles do not match up exactly because one of them does not have a spatial reference...?
(I work with ArcMAP10)

Any hints and ideas what tool will bring me there?
Any help is greatly appreciated
Sandra


texts: How to search for a text from column1 in column2 in QGIS?

texts:
I have street names in a column  (column1) as well as in another one (column2). I would like to write a selection expression to find out if there are similar street names (or parts of the names) between column1 and column2.

I tried the following expression, but nothing happened. Basically, this is what I would like to achieve but it seems like this is not the right solution:

"column1" like '%"column2"%'


texts: Optimizing Intersection Area Loops on big layers using PyQGIS

texts:
Using 3.10.1-A Coruña on Windows 10

I am looping over some (10) national shapefiles and finding intersection and areas but the process is very slow.

I have two loops - one loop for each postcode and one loop for each feature in the Parent polygon layer.

It takes an hour to test 31 Parent polygons against 2,670 postcodes.

It takes 2.25 hours to test 43 Parent polygons against the same number of postcodes.

I have some shapefiles with hundreds and thousands of polygons and I fear this will literally take weeks of computation.

I'm using Ben W's code suggestion from my post PyQGIS layer intersection area seems really inaccurate. Ben W has been a great help to me on this site. I do not want to give the impression that I disrespect Ben W skills and helpful experience at all.

Are there optimization strategies for these kind big loop QGIS tasks?

I'm trying Ben W's suggestion of a spatial index. I've put together this code, but am stuck on the next step:

POA_layer = QgsProject().instance().mapLayersByName('POA_2016_AUST')[0]
parent_region_layer =
QgsProject().instance().mapLayersByName('PHN_boundaries_AUS_May2017_V7')[0]

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

#POA_features = [f for f in POA_layer.getFeatures()]
#parent_region_features = [f for f in parent_region_layer.getFeatures()]

POA_features = {feature.id(): feature for (feature) in
POA_layer.getFeatures()}
parent_region_features = {feature.id(): feature for (feature) in
parent_region_layer.getFeatures()}

POA_index = QgsSpatialIndex()
map(POA_index.insertFeature, POA_features.values())

Parent_index = QgsSpatialIndex()
map(Parent_index.insertFeature, parent_region_features.values())


texts:
I failed to geocode intersection in ArcMap with World Geocode Service provided by ArcGIS Online. The data I have contains street, cross street, city, state. But there is no zipcode

Here is a sample record. The street name here is a combination of street and cross street:

Street: Research Blvd @ Spicewood Springs Blvd

City:   Austin

Region: Texas

Country: US

The weird thing is Geosearch Services can correctly find the address "Research Blvd @ Spicewood Springs Blvd, Austin, Texas, US". But if I try to use the Geocoding Services, it failed. All the addresses in Austin are converted into the same x,y coordinates.


texts: ogrinfo: Get a sample row from a GeoJSON file?

texts:
I know how to examine the properties of a GeoJSON file:

ogrinfo my.geojson OGRGeoJSON -so

But is there a way that I can retrieve a single feature's actual properties, to see what they look like?

I guess I can do this with head, but it would be nice if there was a structured way.


texts:
I have got the following problem:
I would like to access a row from a DBF table. I have tried to access it through SearchCursor, but I am not able to do it without knowing the field name. This is my code:

# Create SearchCursor
rows = arcpy.SearchCursor("tabulka.gdb\zrazkove_uhrny")

# Get names of all fields
field_names = []
sFields=arcpy.ListFields("tabulka.gdb\zrazkove_uhrny")
for field in sFields:
    field_names.append(field.name)

# Print all values of second row
for row in rows:
    print row.field_names[1]

This code works only if I know the name of a row which I don't. I don't even know the number of columns in the table. This code work only if I substitute row.field_names[1] with row.F2 (F2 is the name of the column).

My questions are:

Is there a simpler way how to access values in DBF table?
Is there a way how to access values in DBF table using SearchCursor and without knowing field (column) names?


texts:
I seek for a way to convert, in the print composer, shape into shapefile.
I added rectangle in the print composer with "Add shape" function

I try Rectangles ovals digitizing plugin, but in create shape files outside the print composer. In the ArcGIS, I can do in with "Convert Graphics to features" function


texts: Seeking Python Code Generator in ArcGIS Pro

texts:
I would like to know if there is any way to use Python in ArcGIS Pro in a similar way to SPSS syntax.
I mean if there is a way that while I click in ArcGIS the code starts to auto-generate and then I can copy/edit it the way I want.
I ask this because I know that with R or Python someone can automate processes but I don't know anything about those languages, and maybe this way would be easy to get familiar with the code (like in SPSS).


texts:
I'm not entirely certain, but I think this forum is most appropriate for this question.
I know the location of three observers. A nearby source, who's location is unknown emits a beam of radiation at an unknown time. I know the time at which each observer "saw" the radiation, and I also know the propagation speed of the beam (namely, the speed of light).
Is there a way to determine the location of the source? Or, perhaps, to determine some bits of information about it's location (its direction, etc)? I've been doing some research into multilateration techniques, however I'm unable to find anything that solves this particular problem.
[EDIT]:
This seems to be similar to what I'm after: Multilateration Code
Namely, I want to code an algorithm that, given the timestamps and locations of the observers, can compute the location of many sources. However, the answers to the above question aren't particularly useful.


texts: For loop runs correctly but only on one layer - PyQGIS

texts:
I am trying to automate the creation of null raster files from vectors.
I have this code:
import glob
files = glob.glob("C:/Users/joel/folder/*.shp")
for file in files:
    iface.addVectorLayer(file, "", "ogr")
    processing.runAndLoadResults("saga:rasterize", {'INPUT':file,'FIELD':'max','OUTPUT':'2','MULTIPLE':'4','LINE_TYPE':'0','POLY_TYPE':'1','GRID_TYPE':'3','TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX':None,'TARGET_USER_SIZE':'1','TARGET_USER_FITS':'1','GRID':'TEMPORARY_OUTPUT'})

rasters = QgsProject.instance().mapLayersByName('Rasterized')

for raster in rasters:
    processing.runAndLoadResults("grass7:r.null", {'map':raster,'setnull':'-99999','null':'0','-f':False,'-i':False,'-n':False,'-c':False,'-r':False,'output':'TEMPORARY_OUTPUT','GRASS_REGION_PARAMETER':None,'GRASS_REGION_CELLSIZE_PARAMETER':'0','GRASS_RASTER_FORMAT_OPT':'','GRASS_RASTER_FORMAT_META':''})

I am testing the code on a folder with four shapefiles in. The rasterization process works fine, but when I try to run r.null the process runs four times on the same raster layer when I would like it to run once on each layer.


texts: DynamicCacheLayerManager.Connect deletes my tiles

texts:
You are supposed to be able to connect to your existing tiles with this command.

dynamicCacheLayerManager.Connect(tileCacheDir, tileCachePath);

When I do that my tiles are destroyed and rebuilt.  Clearly this isn't the desired behavior, does anyone have any tips?

I'm using ArcEngine 10 in C# with VS2010.


texts:
Is there a good updated source for comparing existing Geocoding APIs (i.e. Google, Mapquest) and their features?

I'm really interested in their language support as well. For example, Google has libraries for its Geocoding service in JavaScript, Java, Python, and many more.


texts:
we want to use qgis-server within docker behind a reverse proxy.
now qgis-server is running and the get-capabilities works.
the problem is the host, which is given us by the get-capabilities. we get the host and the port from the container and not the host:port which we need for the get map command. i.e.: http://qgis-server.com:8080 is what we have in the get capabilities, but we need: https://qgis-server.com
now there is an option to set the port within qgis-server: https://github.com/qgis/QGIS/blob/master/src/server/qgsfcgiserverrequest.cpp#L48 but this works only if there is not allready a port in the url.
so is there a way to specify the host, port, and protocol of the front facing reverse proxy in a QGIS Server option that QGIS Server generates the appropriate XML?
thanks