texts: How can I set QGIS to not mess up my TIFF images when I geolocate them?

texts:
I am using QGIS 3.10.0-A
I have some TIFF images of simple line maps.  I worked on the TIFFs to reduce file size while retaining quality and got them to where I want (for example, I used Irfanview's Batch Processing with 16 bit/4 colors (black and grays) and LZW compression, 300 dpi, 2947 x 4385 pixels, to get good looking lines with 408 KB file size).  I have tried various ways of geolocating these TIFF files in QGIS and then trying to save the geolocated Tiffs, but every way I try screws with the image and makes it look bad, apparently lowering the number of colors to two (jaggy/pixilated lines like a FAX), or maybe using a different compression scheme (makes the lines look like fractals on the edges).  And the file size goes up dramatically (much lower image quality, much larger file size).
I have tried using the Freehand Raster Georeferencer plugin to do the geolocating, and everything looks great until I use the plugin's "Export raster with world file," which edits my TIFF and makes it look bad (and when I try to open the file in QGIS it opens at 0 lat, 0 long).  If I try to use print composer to save the Freehand Raster Georeferencer layer as a GeoTIFF that also does not work.
I have tried using QGIS's built in Georeferencer to geolocate the TIFFs, which locates them properly, but makes the lines look like fractals on their edges.  I have tried to use print composer to save those rasters as GeoTIFF files, but it also makes the Tiffs look bad.
Trying to find a workaround I have tried to use the World file from the same TIFF, generated by Freehand Raster Georeferencer or print composer, with my good quality Tiff (by making the file names match except for the extension), but the TIFF does not geolocate properly.  I read that if a TIFF has geolocation in its header then the World file is not used, so resaved the TIFF in Paint.net to be sure that was wiped clean if present, but still no joy.
It seems this should not be so difficult, which makes me wonder if there is a setting somewhere for how QGIS handles the TIFF image when it is georeferenced.  To me, the default should be to not mess with the image, as people likely have it like they want it when they go to georeference it.
So, is there a way to set QGIS to simply add the Geolocation information to an existing Tiff file without changing the image itself?
I expect there is a way to do this with GDAL, but I am no wizard and would need clear step by step instructions.


texts: QGIS older version to download onto a Mac ios 10.7.5

texts:
I am unable to download the newest QGIS to my macbook pro (2011, ios 10.7.5) because I need ios 10.8. Unfortunately, I am also unable to update my computer at this time. On th QGIS website there is an option to download older versions, however I don't know which version to use.

Can anyone help/tell me what version I should download?


texts: Using Calculate Field (Data Management) tool to calculate 3D Length in ArcGIS 10.1

texts:
I'm relatively new to the world of model builder and scripting, so it's likely there is a very simple solution to the question I've been trying to figure out.

I'm trying to create a tool that takes a given shapefile, interpolates it over a given surface, adds a "Length_3D" field, then calculates that length (3D) in the attribute table using the Calculate Field tool in the Data Management Toolbox. I've done the same for 2D length using !shape.length@miles! as the expression and I have no issues there.

What I can't seem to figure out is how to write an expression that will allow me to calculate the 3D length from my already interpolated line. Right now the model runs all the way through until this point then throws an error. I'm able to open the newly created ZM shape and calculate the 3D length by using the "Calculate Geometry" tool in the attribute table then selecting "3D Length". I'd really like to be able to perform this in the model though. Has anyone else had a similar issue or know how to write the expression? I've tried "!shape.3DLength@miles!, !shape.3D_Length@miles!, and !shape.3D Length@miles!" with no luck. Thanks!


texts: I have been following the news about the launch of the PAZ satellite (a SAR imaging sensor). The orbit will be sun-synchronous. I know this is important for optical sensors to make sure they have the same illumination conditions. My doubt is why is it relevant for an active sensor? I have heard the same explanation (illumination conditions) but since it is SAR satellite, shouldn't it be independent of that?


texts:
I have code that creates 3 different layers using masking, each time other different pixels get mask and the others get 1 so the final result is 3 layers that are complemntary:

My problem is that when I mosaic those three layers, because each one has only values of 1, when I mosic them I get one layer with the same color and I lose their classes.

My question is if there is any way to determine the unmasked pixels values (for example, instead of 1, that it will be 2,or 3,) and then when I mosaic it i'll not lose the difference colores.

Here is my code, the  part of the masking and the mosaicing:

//create the layers
var imageNDVI=MyImage.select('NDVI');
var gtPOS=MyImage.gt(posBorder).selfMask().rename('range');
var ltNEG=MyImage.lt(negBorder).selfMask().rename('range');
var betMEAN=MyImage.gt(negBorder).and(imageNDVI.lt(posBorder)).selfMask().rename('range');

var ndviClassCol = ee.ImageCollection.fromImages([ltNEG,betMEAN, gtPOS]);

// Mosaic the ImageCollection.
var ndviClassImg = ndviClassCol.mosaic();

// Display the classified mosaic to the Map.
Map.centerObject(geometry, 16);

Map.addLayer(gtPOS,{palette:['#9370DB']});
Map.addLayer(ltNEG,{palette:['fff200']});
Map.addLayer(betMEAN,{palette:['0fff50']});

Map.addLayer(ndviClassImg, {palette: ['ff0000', '00ff00', '006622'], min: 1, max: 3},'mosaic');

My end goal:  to determine the value of the unmasked values

EDIT: I have tried to multiply the values as Daniel suggested, and it worked but the mosaic failed:

This is how I claculted:

//create the layers
var imageNDVI=MyImage.select('NDVI');
var gtPOS=MyImage.gt(posBorder).selfMask().rename('range');
var ltNEG=MyImage.lt(negBorder).selfMask().rename('range');
var betMEAN=MyImage.gt(negBorder).and(imageNDVI.lt(posBorder)).selfMask().rename('range');

var pos=gtPOS.multiply(3);
var neg=ltNEG.multiply(2);

var ndviClassCol = ee.ImageCollection.fromImages([neg,betMEAN, pos]);

// Mosaic the ImageCollection.
var ndviClassImg = ndviClassCol.mosaic();

// Display the classified mosaic to the Map.
Map.centerObject(geometry, 16);

Map.addLayer(pos,{palette:['9370DB']});
Map.addLayer(neg,{palette:['0fff50']});
Map.addLayer(betMEAN,{palette:['fff200']});

Map.addLayer(ndviClassImg, {palette: ['ff0000', '00ff00', '006622'], min: 1, max: 3},'mosaic');

also tried:

//create the layers
var imageNDVI=MyImage.select('NDVI');
var gtPOS=MyImage.gt(posBorder).selfMask().rename('range');
var ltNEG=MyImage.lt(negBorder).selfMask().rename('range');
var betMEAN=MyImage.gt(negBorder).and(imageNDVI.lt(posBorder)).selfMask().rename('range');

var pos=gtPOS.select('range').multiply(3);
var neg=ltNEG.select('range').multiply(2);

var ndviClassCol = ee.ImageCollection.fromImages([neg,betMEAN, pos]);

// Mosaic the ImageCollection.
var ndviClassImg = ndviClassCol.mosaic();

// Display the classified mosaic to the Map.
Map.centerObject(geometry, 16);

Map.addLayer(pos,{palette:['9370DB']});
Map.addLayer(neg,{palette:['0fff50']});
Map.addLayer(betMEAN,{palette:['fff200']});

Map.addLayer(ndviClassImg, {palette: ['ff0000', '00ff00', '006622'], min: 1, max: 3},'mosaic');

as you can see, the layers genereted but the mosaic do not.

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


texts:
When i work on QGIS project (named "pro"), some file with this symbol appears in the workspace folder with the name "pro.qgs~":

I tried to open it- but with no success.

How can i open it and what is the purpose of it?
What to do, so this file will not create in the folder?


texts: (arcengine) export selected features from sde to shapefile (with joined fields)

texts:
I have a SDE featureclass that is joined to an access table in my map control.  I want to export selected features of the SDE layer to a shapefile AND retain the joined fields.
this following method works great (except for retaining the joined fields)

public void ExportMapLyrToShapefile(ESRI.ArcGIS.Carto.ILayer pInLayer, string inputName, string outputPath, string qry)
    {
        // Set the GeoFeatureLayer to the layer in the map so that we can use it as the input feature later
        IGeoFeatureLayer pInGFLayer = (IGeoFeatureLayer)pInLayer;
        IFeatureLayer pInFeatureLayer = (IFeatureLayer)pInGFLayer;
        IFeatureClass inputFeatureClass = pInFeatureLayer.FeatureClass;
        IDataset inputDataset = (IDataset)inputFeatureClass;
        IDatasetName inputDatasetName = (IDatasetName)inputDataset.FullName;

        // Get the layer's selection set (the selected/highlighted features of the layer in the mapcontrol)
        IFeatureSelection featureSelection = (IFeatureSelection)pInFeatureLayer;
        ISelectionSet selectionSet = featureSelection.SelectionSet;

        // Create a feature class name for the output shapefile and open the output's workspace
        IWorkspaceFactory shpWorkspaceFactory = new ShapefileWorkspaceFactoryClass();
        IWorkspace shpWorkspace = shpWorkspaceFactory.OpenFromFile(outputPath, 0);
        IDataset shpWorkspaceDataset = (IDataset)shpWorkspace;
        IWorkspaceName workspaceName = (IWorkspaceName)shpWorkspaceDataset.FullName;
        IFeatureClassName shpFeatureClassName = new FeatureClassNameClass();
        IDatasetName shpDatasetName = (IDatasetName)shpFeatureClassName;
        shpDatasetName.WorkspaceName = workspaceName;
        shpDatasetName.Name = qry.Replace(" = ", "");

        // Use the IFieldChecker interface to make sure all of the field names are valid for the output shapefile.
        IFieldChecker fieldChecker = new FieldCheckerClass();
        IFields shapefileFields = null;
        IEnumFieldError enumFieldError = null;
        fieldChecker.InputWorkspace = inputDataset.Workspace;
        fieldChecker.ValidateWorkspace = shpWorkspace;
        fieldChecker.Validate(inputFeatureClass.Fields, out enumFieldError, out shapefileFields);

        // We also need to retrieve the GeometryDef from the input feature class.
        int shapeFieldPosition = inputFeatureClass.FindField(inputFeatureClass.ShapeFieldName);
        IFields inputFields = inputFeatureClass.Fields;
        IField shapeField = inputFields.get_Field(shapeFieldPosition);
        IGeometryDef geometryDef = shapeField.GeometryDef;

        // Now we create a feature data converter to "Export" the selected features to a shapefile
        IFeatureDataConverter2 featureDataConverter2 = new FeatureDataConverterClass();
        IEnumInvalidObject enumInvalidObject = featureDataConverter2.ConvertFeatureClass(inputDatasetName, null, selectionSet, null, shpFeatureClassName, geometryDef, shapefileFields, "", 1000, 0);
    }

the issue is with the IFieldChecker interface.  it validates all the fields and removed the joined fields before exporting to a shapefile.

Does anyone have any insight on this issue?  If so, please help!  I'm stumpped!

btw... i also tried using an IDisplayTable as my input feature and I get the same output (missing joined fields)

Thanks,
Derek


texts: Overlapping score for each feature in the merged layer after spatial joins

texts:
I've clicked on the Join attributes by location button and can join two layers. But in the merged layer, how do I find out percentage of overlap/intersect calculated for each feature?

For example, in the result, I'd like to know the "overlaps" score for each feature in the merged layer.


texts: getting scratch and current workspace ArcGIS10 using Java

texts: I'm trying to get the scratch and current workspaces returned as path names (Strings) using java in ArcGIS10.  however I can't seem to get anything returned.  I'm guessing it's to do with accessing the environments in Arc but the API is woeful.


texts:
I'm learning how use GeoServer in production. There is a thing that isn't clear for me: where GeoServer saves all the users data?
I know that there is the data folder inside webapps/geoserver/, I use it, but this is the only one place?
For example there is a Disk Quota settings that give the possibility to use PostgreSQL or Oracle instead of H2 DB. Normally I put this option disabled then I don't need to use the H2 DB or other type. How GeoServer remembers the association between user and layers? There is a file or a DB?


texts: Moving GPS Points with a VB script - Bearing works well but distance is wrong

texts:
I am using a script I found online created in 2003 by one of the commentators here.  It is a great script and exactly what I have been looking for and I'm very glad to have found it.  It simply takes a points file and using two fields in the table, bearing and distance it moves the points based on these two fields.  When I run it, the points always seem to move to the correct bearing but the distances are never correct.  For example, I may need to move a point 30 meters away from the original position at a bearing of 85 degrees.  The points move much more than that (hundreds of kilometers or more).  Is this a projection issue?  When I enter 30 in my table as my distance, do some projections interpret this as 30 degrees?  How do I get ArcMap to recognize that this is 30 meters?  If I can solve this I will be extremely happy.  My points file covers very large sections of geography covering many states.  Thank you

The script is below

Option Explicit
Sub MovePoints()
    ' move points in layer(0) based on bearing and distance
    Dim pMxDoc As IMxDocument
    Set pMxDoc = ThisDocument

    Dim pFLayer As IFeatureLayer
    Set pFLayer = pMxDoc.FocusMap.Layer(0)

    Dim pFCur As IFeatureCursor
    Set pFCur = pFLayer.FeatureClass.Search(Nothing, False)

    Dim lBFld As Long, lDFld As Long
    lBFld = pFCur.FindField("Bearing")
    If lBFld = -1 Then
        MsgBox "bearing field not found"
        Exit Sub
    End If

    lDFld = pFCur.FindField("Distance")
    If lDFld = -1 Then
        MsgBox "distance field not found"
        Exit Sub
    End If

    Dim pFeat As IFeature
    Set pFeat = pFCur.NextFeature
    Do Until pFeat Is Nothing
        Dim pConstPoint As IConstructPoint
        Set pConstPoint = New Point
        Dim dAngle As Double, dDist As Double, PI As Double
        PI = Atn(1#) * 4
        dAngle = Azi2Cart(pFeat.Value(lBFld)) * PI / 180
        dDist = pFeat.Value(lDFld)
        pConstPoint.ConstructAngleDistance pFeat.ShapeCopy, _
                                           dAngle, _
                                           dDist
        ' maybe set distance to 0 here
        ' to avoid moving it again
        'pFeat.Value(lDFld) = 0
        Set pFeat.Shape = pConstPoint
        pFeat.Store
        Set pFeat = pFCur.NextFeature
    Loop
    Dim pAV As IActiveView
    Set pAV = pMxDoc.FocusMap
    pAV.Refresh
End Sub

Function Azi2Cart(ByVal dCompassAngle As Double) As Double
    Azi2Cart = -dCompassAngle + 90
End Function


texts: Using ArcMap Field Calculator to count number of identical values?

texts:
I have a file with an attribute table with data in 3 columns, and would like to use field calculator to populate a 4th column. The data is integer numbers between 0 and 56. I want the new column to tell if the numbers in the 3 previous columns are identical or not. I would like the 4th column to state "3" if all 3 columns have the same number, "2" if 2 are the same, and "0" if none of the columns contain the same number.

Does anyone know how to do this calculation?

I tried the following in VB script, but I still get an error message:


texts: No SAGA folder in processing options dialog box

texts:
I am using QGIS 2.18. I am facing this error while using EasyAHP:

  "The error occured during calculation. Please make sure that
  Processing Toolbox's SAGA provider is installed."

But, I guess the correct version of SAGA is installed, but it is not able to recognize the path.

There is no textbox to enter the path.
What should I do now?


texts: Is it possible to map a list with 'set' method for getting directly a list, instead a matrix, avoiding for-loop use?

texts:
The use of for-loops is discouraged in Google Earth Engine because there is pointed out that the same results can be achieved using a map() operation with a specific function independently applied to each element. In following code, I have an initial list of zero values where I want to introduce ones in positions established by a list of indices of non nulls values (idxList). I spent a lot of time trying to find out a map operation for producing that list but, instead a list, it produces a matrix with values placed in adequate positions but it requires additional processing for getting a single list. However, the "discouraged" for-loop produces my desired solution by using directly 'set' method. Functional code looks as follows:
var option = ee.List.repeat(0, 8);

print("initial list", option);

var option2 = option;

var idxList = ee.List([1, 3, 4, 5, 6 ]);

print("index non nulls values", idxList);

var n = idxList.size().getInfo();

for (var i = 0; i < n; i += 1) {

  option = option.set(idxList.get(i), 1);

}

print("desired solution but with for loop", option);

var new_list = ee.List([0, 1, 2, 3, 4]);

var int_lists = new_list.map(function setting (ele){

  return option2.set(idxList.get(ele), 1);

});

print("result mapping a list", int_lists);

Is it possible to map a list with 'set' method (or using its 'iterate' method) for getting directly a desired result list, instead a matrix, avoiding for-loop use?


texts:
I have WMS polygonal layer and want to draw polygon's border with the following style:

My symbol's description:

SYMBOL
  NAME "double_line"
  TYPE vector
  POINTS
    0 0 1 0 -99 -99 0 1 1 1
  END
END

CLASS section in mapfile:

    CLASS
    STYLE
        OUTLINECOLOR 255 0 0
        SYMBOL 'double_line'
        SIZE 3
        WIDTH 1
    END

But in this case you can see self intersection in polygon's edge:

What is the righ way to get correct border style?


texts: ArcPy Feature To Line failing on 'large' datasets?

texts:
I knew that FeatureToLine had issues with large datasets and I have seen other questions revolving around the same error, but I am completely baffled by this. Basically, this function works perfectly with few features (less than a couple hundred), but when a larger number (300+) is entered, it gives this error.

This is the offending function:

    add_print_message("Extracting all intersects as points...")
    arcpy.Intersect_analysis(feature_class, feature_class + "_inter", output_type="POINT")

    add_print_message("Splitting all lines at intersected points...")
    arcpy.SplitLineAtPoint_management(feature_class, feature_class + "_inter", feature_class + "_split")

    add_print_message("Repairing geometry...")
    arcpy.RepairGeometry_management(feature_class + "_split", "DELETE_NULL")

    add_print_message("Converting all split features to lines...")
    arcpy.FeatureToLine_management(feature_class + "_split", feature_class + "_ftl")

It seems to fail after FeatureToLine runs and it gives this error: "Invalid Topology [Feature not found on delete.]"

There are other questions regarding this error on gis.stackexchange but they are of little help.

I have tried Dice and the error still happens and/or Dice fails on chunks that the function utilizes perfectly. I have also tried running it in ArcCatalog but have the same issues.

Any ideas about what could be going wrong (I can't partition any smaller, although I guess I could subdivide, but that is a last resort for now)


texts: Raster clipping over 1300 km2 orthophotos

texts:
I have all of Denmark in 25 cm orthophotos added to my ArcMap (mosaic dataset). I've added a polygon of an area covering over 1300 km2. I've used the raster clip tool to clip the data into a new raster dataset. All data is on a network drive. This process has taken over 32 hours so far.

Is there a difference in speed and such by doing raster clip or is it better to maybe use export data?

Or what is faster? Without moving data locally.


texts: Is there any way for the default text of each annotation subtype to be different? I tried modifying my TextString default value through ArcCatalog but that does not seem to produce any results.


texts: How to convert DHHN2016 to mean sea level?

texts:
I am new to GIS. How can I convert DHHN2016 to mean sea level?

I am using this dataset https://www.opengeodata.nrw.de/produkte/geobasis/hm/dgm1_xyz/dgm1_05314000_Bonn_EPSG4647_XYZ.zip which I stiched together to one geotiff.


texts:
I am trying to cluster points in PostGIS. I have the points and the corresponding Latitudes and Longitudes which I have converted into point geometries. I want to cluster points in such a way that all the points within a cluster are within 25 miles of each other.

My input looks like this:

My desired output is this:

Here the new column Cluster_id denotes the cluster to which the point is assigned to. Please tell me the right query to get this done. I am new to SQL so I am finding hard to write the right query.

So far, I have tried this:

SELECT ST_AsText(unnest(ST_ClusterWithin(the_geom, 1))) FROM all_locations


texts: Transforming point using NetTopologySuite

texts:
I'm exploring how we can use NetTopologySuite to address a geoprocessing task. based on some code sample I found, I was able to piece together a sample application that reads a shapefile and checks if a point falls within a polygon. The issue is, my point is in WGS84 but the shapefile is in EPSG:2810. I need to transform the input point to the CRS of the shapefile.

How do I convert a point from WGS to EPSG:2810 using NetTopologySuite 1.13.3?

Here's some of the code I have so far:

string shpFilename = @"C:\...\polygons.shp";
GeometryFactory factory = new GeometryFactory();
var coord = new Coordinate(-90.00, 45.00);
var point = factory.CreatePoint(coord);

I should also note that when I manually convert from WGS84 to the projected system and use those converted values, my code works as designed. This confirms this is a matter of transforming the coordinates before attempting to determine if they fall within a polygon.


texts: GRASS GIS using a loop with r.neighbors command

texts:
I'd like to do a trend analysis with ndvi data from 1982 til 2005. I am using the r.neighbors command to get first the spatial variance (then I would use r.series method=slope to get the actual trend) and in order to speed things up I'd like to do it using a loop... however it is not running :(
Here is my code... can someone help me?

for a in
  1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993
  1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 ; do

  for i in 311 ; do
    r.neighbors
      input=AVHRR_NDVI_$aDOY$i
      output=AVHRR_NDVI_$aDOY$i.sp.var
      method=variance
      size=5
     echo "loop ends" ;
   done ;
 done


texts:
Does an INSEE code exist for French regions?
Is there a list of the official INSEE codes for French regions?


texts:
In a MapBox Studio project of mine, I have field values that are stored as JSON values.  My features all have attributes like:

id 48545756
meta {}
relations []
tags {"access":"forestry","highway":"track","source":"WRI;DIAF","start_date":"before2013-03-28"}
type way

I'm trying to filter by a key in the tags field.  Something like:

#layer['tags.access' = 'forestry']{ ... }
or
#layer['tags' = '{"access": "forestry"}']{ ... }

But of course MapBox Studio interprets the JSON as a string, not something that can be parsed.

Any thoughts on how to get around this?

Why JSON as feature attribute values?  MapBox's uploads api provides a really good way to programmatically upload data to a MapBox account.  If you send a geojson file (rather than a mbtiles file), MapBox tiles it out with it's attributes as a json value in the 'tags' field.


texts: How can I get a QGIS2Web project to display info for multiple layers in a popup?

texts:
I have built a project in QGIS and QGIS2 Web that has multiple polygon layers that overlap each other. I am trying to get the info pop up to display the data for all of the overlapping layers at a point.

I have tried exporting the project to OpenLayers and to Leaflet neither seems to work.

The OpenLayers version is here http://knowwhereconsulting.co.uk/WebMaps/lgsstest/
The Leaflet version is here http://knowwhereconsulting.co.uk/WebMaps/LGSS/

Anyone know how to fix this?


texts:
I have a grid on a map of my city. I have 1500 plant species in a CSV and need to print each species on a copy of the map in the grid. So I will have 1500 PDFs of the same map but with a different point layer representing which grid squares contain each species.

Is there a way to automate this somewhat so I don't have to create 1500 different print layouts? I can have each species as a layer or names as attributes in 1 layer, there are only:X,Y coords, species name data and red list status as data.


texts: Stratified Random Sampling by ID in QGIS?

texts: I've been playing around with the Random Points tool in QGIS to realize stratified random sampling, but I'm running into a problem where I want my strata to be based on a polygon ID, not each individual polygon.  In other words: I may have multiple non-contiguous polygons that all represent a single strata, where I want a TOTAL of (e.g.) 10 points randomly distributed among the disconnected polygons, not 10 per polygon.  Is there any way to do this?  Would a multi-part polygon work, for this?


texts: QGis 3.16 Script exported from MODEL cannot save results

texts:
I built a simple model to perform CLIP operations, with 3 inputs: Input features, Clip features and destination folder determined by the user, as shown in the image below:

In "Output features" of File/Folder input is written "@TestFolder" (which is the path to the user selected folder), plus the resulting file name ("clipped.shp").
Everything works fine, and the result shape is saved to the path chosen by the user.
However, when this model is exported as a python script, no result shapefile is saved to the chosen folder; bellow is a sample of this python script:
class ModelTest(QgsProcessingAlgorithm):
def initAlgorithm(self, config=None):
    self.addParameter(QgsProcessingParameterVectorLayer('Inputpolygon', 'Input polygon', types=[QgsProcessing.TypeVectorPolygon], defaultValue=None))
    self.addParameter(QgsProcessingParameterVectorLayer('Clippolygon', 'Clip polygon', types=[QgsProcessing.TypeVectorPolygon], defaultValue=None))
    self.addParameter(QgsProcessingParameterFile('TestFolder', 'TestFolder', behavior=QgsProcessingParameterFile.Folder, fileFilter='All files (*.*)', defaultValue=None))

def processAlgorithm(self, parameters, context, model_feedback):
    # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
    # overall progress through the model
    feedback = QgsProcessingMultiStepFeedback(1, model_feedback)
    results = {}
    outputs = {}

    # Polygon clipping
    alg_params = {
        'CLIP': parameters['Clippolygon'],
        'S_INPUT': parameters['Inputpolygon'],
        'S_OUTPUT': QgsExpression(' @TestFolder + \'/ClippedPolygon\'').evaluate(),
        'S_OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
    }
    outputs['PolygonClipping'] = processing.run('saga:polygonclipping', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
    return results

It is something missing here, for sure, so I ask for some help.


texts: Off by about 85 meters when converting NAD27 UTM to WGS84 lat/lon

texts:
Looking through other answers on converting to or from NAD27 CONUS, I don't
see an answer that seems to fit my problem.

I am creating paper topographic maps by layering information derived from USGS DEM data (NAD83) for contours and slope shading. I am using OpenStreetMap (OSM) data (WGS84) for trails, roads, water features, etc. Finally I am overlaying a UTM grid. All this is being formatted using Mapnik.

Everything is aligning well (trails and road appear correct with respect to the contours. The waterways from OSM also nicely match up with the drainages as shown by contours derived from USGS DEM data, etc.) And if I use WGS84 for my UTM grid it seems to line up well too, at least the positions of map features are the same within my measuring capability as a USGS topo map. But if I overlay a NAD27 UTM grid my easting grid lines appear to be off by about 85 meters. The northing grid lines appear to be correct. I could live with 5 or 10 meters as that is the accuracy of my GPS but not 85 meters.

I create the UTM grid by determining the UTM measurements for the edges of my map and then going to points the nearest even 1000 meters off paper and putting ways/lines between the points so located. The position of the points is converted to from UTM to lat,lon and the points and ways are written to a OSM XML file which is used as the source for the UTM grid layer in Mapnik.

Conversion of UTM to lat,lon is using pyproj in a Python script. My instantiation of the projection objects is as follows:

wgs84utm = pyproj.Proj(proj='utm', zone=v['zone'], datum='WGS84', ellps='WGS84')
nad27utm = pyproj.Proj(proj='utm', zone=v['zone'], datum='NAD27', ellps='clrk66')

When I convert a point from UTM to lat,lon I use the following:

if (proj == 'wgs84'):
    point = wgs84utm(easting, northing, inverse=True)
elif (proj == 'nad27'):
    point = nad27utm(easting, northing, inverse=True)

One thing that confuses me is I only see one NAD27 UTM setup in the pyproj espg list. Based on my GPS and other mapping programs I was expecting to see several (CONUS, Alaska, Canada, etc.). There is a NAD27 CONUS but it is an Albers projection, not a UTM projection. Based on values for that I tried

nad27utm = pyproj.Proj(proj='utm', zone=v['zone'], datum='NAD27', ellps='clrk66', lat_1=29.5, lat_2=45.5, lat_0=23, lon_0=-96, x_0=0, y_0=0)

But the results were the identical.

I am a newbie with python and map projection conversions, so I'm pretty sure that I've made some mistake in my usage of the projection library but don't know what. I'd appreciate any suggestions on where I have gone wrong.
Thanks!


texts: Remove "Save As" from QGIS Context Menu Using PyQGIS

texts:
I want to remove an existing option ("Save As") from context menu of QGIS layer, what I have found so far is, to add new action in the LegendItems, or overwrite completely with new options. But could not find any option using that I would be able to remove an existing default Context menu option.

What I have also seen that there is menuProvider() method available, but could not get any details with this method.

How can I remove existing "Save As" option from the context menu of QGIS?


texts: Using If Statements in ArcGIS Calculate Field?

texts:
I've never used Python scripting before and can't seem to figure out how to use an If statement in the Field Calculator tool.

I have a field called 'MAX_GRIDCODE' which returns the highest force a piece of infrastructure is exposed to. I need to use this number to calculate a simple estimate of damage in the field 'DI'.

The statements I need are:

If "Max_GRIDCODE" >0 and <=0.5 then "DI" = 1
If "Max_GRIDCODE" >0.5 and <=1 then "DI" = 2
If "Max_GRIDCODE" >1 and <=1.5 then "DI" = 3
If "Max_GRIDCODE" >1.5 and <=2 then "DI" = 4
If "Max_GRIDCODE" >2 then "DI" = 5

Does anyone have tips on how to code this using ArcGIS 10.2 for Desktop?


texts: Disable showing content while moving features in QGIS

texts:
When moving features in edit mode, normally it shows features content (new position of features) in every displacement of mouse cursor as shown in the image.

In case of small number of features, no problem. But while moving thousands of features, QGIS freezes for a long time before it shows features content. Sometimes I have 10.000-50.000 features (~1.000.000 features in total) that I have to move at the same time.

I need to move features manually. I don't move all features at once, but partially. Using tools slows me down.

How can I disable showing content while moving features?


texts: Transparent custom legend with mapbox/leaflet.js

texts:
I am building a map with mapbox.js.  I have created a custom legend and want to make it partially transparent, but can't figure out the right CSS.

Here is my legend html snippet:

My title

My sub-title

I can't figure out the right CSS to make it transparent.

.legend .legend-title {
    text-align: Center;
    margin-bottom: 8px;
    font-weight: bold;
    font-size: 120%;
    color: #ffffff;
    }

.legend{
    background: rgba(80,80,80,0.3);
       }

I have also tried to change the style of the map-legend directly.  I don't think that this is a good way to go at it, but I found this obj in the browser dev tools.

.map-legend.wax-legend {
    background: rgba(80,80,80,.1);
    }

Can anyone help me get this working?

My preference is to use the alpha in rgba on the background and leave the text fully opaque.


texts: What are OpenStreetMap license requirements?

texts:
We are planning to mix OpenStreetMap (OSM) layers with our own layers and afterwards to build an application which will perform common GIS tasks.

Assuming we are allowed to mix our layer and OSM layers in one application, where is the most appropriate place to mention OSM copyright info?

Perhaps including this info in the layer metadata is enough or maybe we have to include this info in the application as well.


texts: How to restrict the creation of vectors (end vertices) to existing vertices in QGIS

texts:
My specific interest lies in network creation, however, knowing how to do that in general would be a plus.
Is it possible to prevent a user from creating a line (line vector) without connecting it perfectly to nodes (point vector) on both ends?

This would save tremendous time in topology checks and corrections later.

I know about using snapping, but it is sadly not fool proof.


texts: SLD file not being read from ArcGIS Server WMS service

texts:
I am trying to link an external SLD file to my WMS service in ArcGIS Server but it is being ignored.
I have one layer in my MXD called 'Rivers', I then published the layer as WMS service in the ArcGIS Server via ArcMap.
I have checked 'WMS' in the capability tab. In ArcGIS Sever Manager, I navigate to my Service--> Capabilities --> WMS --> SLD Path or URL. I have entered my SLD file path, then restart the service. But, change is not reflecting in the service, it is still using ArcMap style.
I have gone through this link, but no success.
My service name is "RiverService". I am using version Arc V. 10.8.
Here is my SLD file:
<sld:StyledLayerDescriptor xmlns="http://www.opengis.net/ogc" xmlns:sld="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0"
    xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
    <sld:NamedLayer>
        <sld:Name>Rivers</sld:Name>
        <sld:UserStyle><sld:Name>Rivers</sld:Name>
            <sld:Title>lineSymbolizer</sld:Title>
            <sld:FeatureTypeStyle>
                <sld:Rule>
                    <sld:LineSymbolizer>
                        <sld:Stroke>
                            <sld:CssParameter name="stroke">#FF00FF</sld:CssParameter>
                            <sld:CssParameter name="stroke-opacity">1</sld:CssParameter>
                            <sld:CssParameter name="stroke-width">2</sld:CssParameter>
                        </sld:Stroke>
                    </sld:LineSymbolizer>
                </sld:Rule>
            </sld:FeatureTypeStyle></sld:UserStyle>
    </sld:NamedLayer>
</sld:StyledLayerDescriptor>


texts: Field calculator uses wrong fields when working with joined data

texts:
I am having an issue in qgis table joins.
i have a joined a csv file to a sqlite database and common id is used for join.
the join is displays correctly in  attribute table.

the attributes in sqlite file is id and the csv file has id, building, surveyorid, surveydate. I want to add the surveydate information to the sqlite data permanently.

In the field calculator. i created a new field survey_date.

the expression field is "surveydate" (field from csv file), the output preview displays the result correctly.

The problem is on clicking ok, the new field(survey_date) gets populated with values of surveyorid. I tried other cases and each time, the new field gets populated with values from the field that is before the field in expression, i.e. if i try with surveyorid, i get the building values.


texts: How does Geocortex Essentials handle projections?

texts: In ArcGIS Online, the projection is defined as whatever the basemap is... Which is typically web mercator auxillary sphere. In Geocortex Essentials, is the projection defined the same way (whatever the basemap is?) and in terms of best practice, should the projection be web mercator auxillary sphere since it is a web mapping application? I would prefer it to be in NAD83


texts: Does QGIS georeferencing set to magnetic or grid north for British National Grid?

texts:
I have just been georeferencing a map on QGIS using British National Grid and it is not lining up with the vector maps I downloaded from OS. I used the 'From map canvas' option and used 'Linear' transformation. However, there appears to be a twist of a few degrees with the map and this suggests that the one is setting as grid north and the other magnetic which would explain the difference. I tried a few times with different reference points, but always the same result.

Anyone any ideas?

Ok, just to add. I rotated the original image by 1.5 degrees and then re-georeferenced and it near enough lined up. Surely the idea of georeferencing would be to put the points in the right place no matter what the orientation of the original raster map?


texts: Why is quantile reclassification method yielding classes of vastly different cell counts?

texts:
I am trying to reclassify many euclidean distance layers (which are based on lines) each into 10 classes of equal area. When I attempt to do this, using the quantile reclassification method, some classes contain twice as many cells as others. I thought that using the quantile method, the cell counts of each class were supposed to be roughly equal. The "Slice" tool does  generate classes of equal cell counts when using the "equal area" method, but I need to know the break values, which the slice tool does not give (But reclassify does). How can I reclassify these layers into ten classes of equal area and know the break values of each class?

I have found that this can be done by first converting the floating point raster to an integer using the floating point tool.


texts: Reprojecting: trouble converting MapInfo data to shapefiles in QGIS

texts:
I was given a dataset generated in MapInfo. I need to use these files in QGIS, so I imported the .TAB files and saved them as ESRI shapefiles. I then realised the original dataset uses a custom "user defined" CRS. When I change the CRS from this user defined system to WGS84 in QGIS, serious positional errors arise! The shape is also wildly altered.

I tried converting the MapInfo data to ESRI shapefiles using ogr2gui and got the same results.

How can I reproject the data to overcome the positional and shape errors? Is there a way of changing the original MapInfo CRS to WGS84, ensuring the shape and desired positioning will be correct in QGIS?


texts: Unary union/Dissolve returns null geometry in QGIS 3

texts:
I am writing a python processing script in QGIS 3. In it, I have polygon geometry objects stored in a list polygonGeoms and I want to dissolve all of it into one geometry. However, when calling the function QgsGeometry.unaryUnion(polygonGeoms), it returns a null geometry. I've used the unaryUnion function before so I know that it should give me the result I want.

To investigate, I exported the geometries as a temporary shapefile layer. And when I use the Dissolve tool in QGIS 3, it also returns a null feature. The Dissolve tool works on the shapefile when I use QGIS 2.18 though, for some reason.

This is a link to the sample shapefile.


texts:
I have a simple model that calculates the NDVI from bands 3 and 4 of a raster dataset.  The input bands are 8-bit unsigned, whereas the default output NDVI is signed 32 bit.

How can I specify the pixel depth for the output NDVI as 8 bit unsigned?


texts:
me and my collegue are making a web map application for my town for tourism purpose. The map itself is pretty done, actually, so we were wondering whether add some photos of the town, or that shows photos of the feature that has been selected. Something pretty much like Google Maps.

Is there any plugin that fits correctly with OpenLayers that does this job correctly?


texts: How to compile Degrib on Ubuntu 12.04 LTS?

texts:
I'm trying to compile degrib on a typical Ubuntu 12.04 LTS Server. However, I get an error "undefined reference to gzopen64" which appears to be common when compiling degrib. Do you have any suggestions?

I see Salda got it working at How to convert GRIB-1 to GRIB-2?. Any suggestions will be much appreciated.


texts: Slicing all values in field using arcpy.da.UpdateCursor

texts:
I have a folder with a (3) shapefiles. Each shp. has a field titled "country" with a country name. I've already set my desired workspace correctly and everything else in my script works for what I want to do, it's just this part I'm struggling with. What I'm attempting to do is to slice the country name so only the first letter of the country is displayed in the country field for each attribute. So for example, if my first attribute in the country field is Canada, I want to see 'C' instead, if my second row is Portugal, I want to see 'P' etc. I need the loop to go though each shapefile in my folder as you can see. I generally don't have much experience at all using string slicing or string concatenation. My desire is to use an update cursor to accomplish this.

Below is a snippet of my current script attempting to accomplish this.

   try:
       for fc in arcpy.ListFeatureClass():
           with arcpy.da.UpdateCursor(fc,["country"]) as cursor:
               for row in cursor:
                   row[0] = [0]
                   cursor.updateRow(row)
               del row, cursor
   except:
       print ("Could not complete the string slicing notation.")


texts: Avoid overlap in QGIS sometimes works and sometimes doesn't for no apparent reason

texts:
I have it set to "avoid overlap on active level" and "self-snapping" and "snapping on intersection" are both selected as well. I read this question: Avoid Overlap option not working in QGIS but it didn't seem to arrive at an answer. I'm very, very new to this and basically all I know how to do is make polygons. As far as I can figure there is no difference between the ones that successfully avoid overlap and those that don't. It happens within the same shp file, the same work session. I draw one polygon and it doesn't overlap, draw another one right next to it and it does overlap.
I've also read this: https://issues.qgis.org/issues/2921 but I really don't understand what it's talking about.


texts: Parsing entire ArcSDE database searching for field for export?

texts: I'm wanting to parse through an entire ArcSDE database (feature dataset, features classes, tables, everything...) searching for a field (compkey) and export the numbers to a txt file.


texts: Getting Uncompressed Size of Raster using ArcObjects?

texts: How can we get the Uncompressed Size of a raster as shown in the layer properties below using ArcObjects?


texts: Convertion of spectral indices formulas from Landsat to Sentinel

texts:
There are many spectral indices in literature. However, the formulas mostly refer to landsat imagery, and thus the corresponding spectral bands.

Question:
How can a formula that was originally designed for Landsat be converted into a Sentinel formula?

Problem:
Not all Landsat and Sentinel bands seem to correlate well, as Sentinel has a higher number of bands that often also differ in the spread of wavelengths covered. These differences have to have some influence on the results when adjusting the formula from Landsat to Sentinel. I'm especially confused in Landsat8 Band 5, which is covered by Sentinel Band 8 and 8a. Furthermore, Sentinel data has Bands 5,6 and 7, which are not covered in Landsat8 data. Older Indices are missing these spectral bands completely.

In some bands they also strongly differ reflectance.

https://scientiaplusconscientia.wordpress.com/2017/03/04/remote-sensing-comparison-landsat-sentinel-visible-infrared-spectrum/

Example:
Let's take the Enhanced Built-Up and Bareness Index (EBBI) formula as an example for discussion (As-syakur 2012).
Band 4 and 6 are easy to pick (Sentinel 4 and 11, respectively)

Next, the formula states to take the Landsat 5 Band. When converting it to Sentinel, the question arises whether to take Band 8 or Band 8a.


texts: Python toolbox field value precision - selection based on Double field

texts:
I have a python toolbox (.pyt) where the user selects a layer, a field from that layer and then values from that field. The tool uses this to make a selection based on the field values. It's working fine for all other datatypes (integer, text/string, float) but there is an issue with double that I think I've narrowed down to precision.

Manually making the selection by attribute the field value is 0.1000000014901162 (works). The value displayed in the tool is 0.100000001490116 (see image).

The field is data type double with a precision of 0 and scale of 0 (see image).

I would like a solution that will work across any precision/scale. I've tried using python's decimal, but that gives too high a precision 0.100000001490116119384765. This one selection is part of a very long script but I've disassembled the functions and here are the relevant bits (not indented/functional):

    #def getParameterInfor(self):
    # FC
    landUse = arcpy.Parameter("Landuse/Greenspace Polygons",
                              "land_use", "GpFeatureLayer", "Optional", "Input")
    # Field
    LULC_field = arcpy.Parameter("Greenspace Field", "LULCFld", "Field",
                                 "Optional", "Input")
    # List of fields from table
    landVal = arcpy.Parameter("Greenspace Field Values", "grn_field_val",
                              "GPString", "Optional", "Input", multiValue=True)

    LULC_field.parameterDependencies = [landUse.name]
    landVal.parameterDependencies = [LULC_field.name]
    landVal.filter.type = 'ValueList'

    #def updateParameters(self, p):
    # Conditions must be met for each input to become active
    if True in [p[4].value, p[6].value]:
        p[17].enabled = True
    else:
        p[17].enabled = False
    if p[17].altered:
        p[18].enabled = True
    if p[18].altered:
        in_poly = p[17].valueAsText
        TypeField = p[18].valueAsText
        p[19].enabled = True
        p[19].filter.list = unique_values(in_poly, TypeField)

    #def main()
    landuse = params[17].valueAsText
    field = params[18].valueAsText
    fieldLst = params[19].values

    # Convert the list of unicode to float
    fieldLst = map(float, fieldLst)

    # Make selection
    arcpy.MakeFeatureLayer_management(landuse, "lyr")
    # Create where clause from list
    exp =''
    for item in fieldLst:
         if type(item) == float:
              exp += '"{}" = {} OR '.format(field, float(item))
    whereClause = exp[:-4]
    sel = "NEW_SELECTION"
    arcpy.SelectLayerByAttribute_management("lyr", sel, whereClause)

I haven't included all the code to differentiate datatypes, but the above currently works on float and other datatypes are differentiated either when converting from unicode or when creating a where clause from the list.


texts: How the create the nearest point on a line layer in QGIS?

texts:
How can I generate a projection of points on the nearest line?
I have a layer with several points and I have another layer with several lines.
And I want to generate a projection of these points but creating points in these projections.
I found this picture online that illustrates what I want to do:

I want to create this projection but also create a point in the final projection (in the line) as a new layer with these points.


texts:
I want to create a buffer of 5km around the coastline, and save that as a separate shapefile.

I have a shapefile with almost 180,000 features. It's a line shapefile, representing the coastlines of the world (I got the data here). The shapefile is in CRS (EPSG:4326, WG84).

So far I have been testing three approaches and all are either very slow or crash:

1) QGIS 2.14.2 Essen (with functioning GRASS functions)
I go to 'Vector' --> 'Geoprocessing Tools' --> 'Buffer(s)...'
I select the shapefile,
For buffer distance I use: 0.05 (EPSG:4326 which is about 5.5km)
I check the 'Dissolve buffer results' options.
This option doesn't seem to work, as QGIS crashes. I expect it has to do with a shortage of memory available. "Error: minidump written to (filepath)"
[

2) PgAdminIII 1.20.0 with the following state I try to create a new table and fill it with the buffer:

CREATE TABLE public.lines_buffer AS SELECT ST_Buffer(geom,0.05) FROM public.lines

But the connection keeps on getting lost and thus breaking the operation. The furthest I got was 2 hours of query time, but that wasn't enough.

3) Python 2.7 (shapely, fiona, descartes) My third attempt, is still in progress. I am trying to loop through the large shapefile and make a buffer. I'm not sure if my code is very efficient. I still have to "dissolve" after the buffer, which would also take a while. I plan on testing shapely.ops.cascaded_union.
EDIT: Even running the first 2 items takes very long...[Spyder IDE: 100% memory usage]
My current loop just limits the first 3 items in the fiona.collection
Credits go to Tom Macwright

def readshape(pathfolder, filename):
filepath = os.path.join(pathfolder, filename)
with collection(filepath) as Input:
    schema = Input.schema
    with collection('buffer_shape.shp', 'w', 'ESRI Shapefile', schema) as output:
        for i in range(0, 3):
            output.write({
                'properties': {
                    'source': Input[i]['properties']['source']
                },
                'geometry': mapping(shape(Input[i]['geometry']).buffer(0.05))
            })

Does anyone have any tips,  which software is the quickest and most stable option?


texts: Output formats not showing in heatmap plugin (QGis 2.10.1 in ubuntu 14.04)

texts:
I just upgraded to QGIS 2.10 (using the debian repos) and the heatmap plugin doesn't show any formats in the dropdown menu.

I'm able to open every raster format I have at hand and, of course, the drivers are listed in the gdal driver options.

I looked around and found a similar question here but refering to QGIS 2.8 and without any answers.


texts: Web app search field with 'meaningful' OSM places

texts:
I'm building a web app based on OSM places, excluding administrative zones (such as countries, states, and counties).  I've taken an extract of place type = city or town or village or hamlet.  But I've noticed a lot of places I need are not there.  However when I search them on OSM's website, they exist.  Example 1: Westport, Massachusetts.  This is a town and has a boundary in OSM but it would not exist when extracting place=town from OSM.  Example 2: Nantucket.  This is a little more abstract... it is an island but is still something people think of as a place so they would search for it.

I've noticed web sites like FourSquare that use OSM have solved this problem.  When you search on FS, the first field is what you are searching (eg Joe's Pizze) the second is place (this is the OSM places).  For place they figured out how to translate OSM places into places someone would search for.  For example I find all the cities/towns I already have but also Nantucket and Westport, MA.  How do sites like FS do it?  Is there an logical approach?

Just like FS, I would want an autocomplete field that a user would search these places through.


texts: Lining up polygon legend patches in ArcGIS for Desktop?

texts:
I have tried the properties and everything, I can't understand why my polygons are not lining up correctly in the legend.

I created this legend the usual way by adding my three layers to the legend wizard and creating the legend this way. It is only one legend, I've just recreated it and the same thing happens.

This is how I am creating my legend:

Insert>legend>legend wizard (select 1 layer)> click defaults


texts: Using Python to send arguments to LAStools basics

texts:
I'm currently trying to develop a LAStools specific routine, setting all the parameters for lasheight then using the output of that to run lasclassify. The closest thing I've found to being able to write this in Python is in Sending arguments to LAStools function from Python script not working

In the above post I'm a bit curious where the input and output .las files are. Also, not exactly sure how that Bat functionality is working. Does anyone have an example of a complete script they would be willing to post?

Here's what I have so far, I can't seem to get it to work though. I think it has something to do with the environments setting. I'm running the exact text from the bat file in the command line and it's working now that I updated the default PATH to where my LAStools bin folder is.

import os, random, sys, subprocess

input_path = "D:\Working_D\LAStools_scratch\input\\*.laz"
height_path = "D:\Working_D\LAStools_scratch\height"
arg_cores = "-cores 16"

directory = r'D:\Working_D\scripts\BatchFiles'
with open(os.path.join(directory, 'BatchTemp.bat'), 'w') as OPATH:
      OPATH.writelines(["lasheight ",
                        "-i {0} ".format('"' + input_path + '"'),
                        "{0} ".format(arg_cores),
                        "{0} ".format("-odir"),
                        "{0} ".format('"' + height_path + '"'),
                        "{0} ".format("-olaz"),
                        "\n "])

BatFile = 'D:\Working_D\scripts\BatchFiles\BatchTemp.bat'

proc=subprocess.Popen(BatFile,stdout=subprocess.PIPE,stderr=subprocess.PIPE,shell=True)
output,error=proc.communicate()

os.remove(BatFile)

print 'process complete.'

SOLVED: I was running this in IDLE and it would not work but running the Python script through the cmd window (e.g. D:\Working_D\scripts\lasheight_v2.py) and it worked!


texts: Margins and adding data to layout in ArcGIS Pro?

texts:
Just starting out with Pro trying to put together a simple map. Is there a way to add printer margins in a layout?

Also I'm confused on how to import my data into the layout itself. I have it in another active window.


texts: Setting up a dynamic table Excel to ArcMap

texts:
I am trying to create a GIS map that contains data for over 200 monitoring wells that are sampled 4 times a year. I would like to set up the GIS map so the attributes table for the well layer updates automatically when I update the excel spreadsheet that I would like to link to the well layer attributes.  Is there an easy way to link an excel spread sheet to an attribute table so the attributes will update automatically when the excel file is updated. I found some information about this topic in the question pool, but didn’t find a clear solution.

I'm using Arcmap by Esri version 10.3?


texts: I am using PostGIS GUI to import shapefile. The connection is good. While I successfully imported the Alabama table. I got the following error when I tried to import the other two shapefiles. I imported those files separately. Why is that?


texts: Why does arcpy.da.SearchCursor() class not autocomplete within Eclipse?

texts:
I am very new to GIS.

I just installed eclipse luna with pydev plug-in and ArcGIS 10.1 for Desktop.

When I write a script, trying to import arcpy.da.SearchCursor I see that no SearchCursor is available to autocomplete.

What is the problem? What I did wrong?


texts: Making union of some raster coverages without the area between them

texts:
I have a set of raster coverages. Some of them are really far from each other. When I use ST_Union function to union them, the function fill the area between them.
I explain it in a simple example:
I have a grid of polygons and select 2 of them (by their id) that are not neighbours.

Now, I have two different polygons in a table. Then, I used ST_AsRaster() to make two rasters out of the two polygons.
At the end I used ST_Union() to make a union of those rasters.

select *
INTO test1
from grid
where id= 277765 or id = 277762;

select st_asraster(tile_border,100,100, '8BSI'::text, 0, -1) as rast
into test1_rast_1
from test1 ;

select ST_Union(rast)
into test1_rast_1_union
from test1_rast_1;

I want the area between to be no-value but the ST_union fills it with 0. The created raster is as follows:

My real data consists of the raster coverages that are really far from each other and I have many layers that I want to union them.
How can I use ST_Union in a way that it does not fill the area between the two geometries? because I think filling it for each iteration takes a long long time.


texts: Uncaught TypeError: o.intersection is not a function in turf.js

texts:
I'm working on a buffer analysis and then intersection using turf.js. I have a buffered output layer to which I need to intersect another polygon layer and get the output. Below is my code,

        function calculateBuffer() {
            var fc1 = turf.featureCollection(fc1Layer.getLayers().map(function(f){
                    return f.toGeoJSON()
                }));

            var buffer = turf.buffer(fc1, 5, 'kilometers');
            buffer.properties = {
                "fill": "#6BC65F",
                "stroke": "#25561F",
                "stroke-width": 2
            };
            bufferLayer.setGeoJSON(buffer);
           //var bufferGeojson=bufferLayer.setGeoJSON(buffer);
    //polygon layer//
            var fc = turf.featureCollection(fcLayer.getLayers().map(function(f){
                    return f.toGeoJSON()
                }));
    //buffer layer//
             var buffer = turf.featureCollection(bufferLayer.getLayers().map(function(f){
                    return f.toGeoJSON()
                }));
  //intersection//
            var polygon = turf.intersect(fc,buffer);
        polygon.properties = {
            "fill": "#25561F",
            "stroke": "#25561F",
            "stroke-width": 10
        };
        intersectLayer.setGeoJSON(polygon);
        }

I'm getting the following error:turf.min.js:16 Uncaught TypeError: o.intersection is not a function. Could anyone help me on this?


texts: I'm using the most up to date version of QGIS. I have created a vector layer of multiple polygons that form areas of one larger area. I cannot get the polygons to connect at the corners or the edges.


texts: Convert netcdf4 to ASCII using Python and gdal: No output

texts:
I want to convert a netCDF4 dataset (including grouped datasets) to Arc/Info ASCII Grids, using Python's gdal package. I tried to orientate at this question, the only thing I added is a loop over the subdatasets. Now, however, something went wrong and the script doesn't produce any output.

Can anyone help? Here's my code:

#Open existing dataset
src_ds_list = glob("..\\*.nc")

for path in src_ds_list:
    src_ds = gdal.Open(path)
    subdata = src_ds.GetSubDatasets()

    for s in range(len(subdata)):
        sds = gdal.Open(subdata[s][0])

        dst_path = '..\\' + subdata[s][0] + 'arbitrary' + '.asc'

        # Open output format driver, see gdal_translate --formats for list
        oformat = "AAIGrid"
        driver = gdal.GetDriverByName(oformat)
        # Output to new format
        dst_ds = driver.CreateCopy(dst_path, sds, 0)
        dst_ds = None
        sds= None
    src_ds = None


texts: iterate over cluster in markerCluster : leaflet

texts:
i need to iterate over all clusters that are currently present in markercluster using leaflet. there's a method called markercluster.getClusters() in googlemaps to do this job.

is there something similar in leaflet too?


texts: Is it possible to automate sewer service line creation?

texts:
I am a GIS Tech working for a utility company, part of the asset maintenance program entails running a camera down gravity mains (sewer pipes) from the upstream manhole to the downstream manhole.

Every time the camera finds a sewer lateral (service line) it records its distance and the angle it leaves the gravity main.

I have been able to use ArcMap to plot the distances as points by converting the gravity mains to a route layer and using the create route events too on the original excel file given to me which contains the distances for each service line.

Now here is my main question - most of our gravity mains currently do not have the sewer laterals (service lines) digitized yet. With the location of the sewer lateral (point file) i also have the direction it leaves the gravity main.

Is it possible some how in arcpy to automate the creation of the sewer laterals, for example, just make all the sewer lines 20 ft and set a predetermined direction (the original points have an angle measurement for the sewer lateral direction leaving the gravity main).

Could anyone give me any direction on how I might go about automating the creation of the sewer laterals


texts:
I am testing Lizmap on a project where I display a high density of points from a point layer.

Anyone knows a way to cluster the points in lizmap or qgis server?

The idea is to end up with some clusters similar to this:

https://www.mapbox.com/mapbox.js/example/v1.0.0/leaflet-markercluster/
http://trac.osgeo.org/mapserver/attachment/ticket/3700/cluster.png


texts:
Is it possible to use Python 2.6 (or 2.7?) with ArcGIS 9.3?

The official stance is no, not supported. There are a number of reports in the wild though that it can be done. How is this accomplished? and is there anything to watch out for?


texts: Wrong vector selected mouse-clicking on computed layer

texts:
[qGIS 3.10.14]
Only on computed layer such as virtual or spatialite view:
When simply clicking in canvas (or using any other gui/spatial selection method) to select particular feature(s), what gets highlighted/selected are never the feature(s) actually at the cursor position, but features situated elsewhere.
It is maddening. Right now I'm just trying to select and copy to a working layer one property parcel from a spatialite view that links 60+k polygons with appraisal attributes. But clicking on the parcel highlights and selects in the attribute table a parcel that sits clear across the county. Clicking and dragging or selecting by polygon/freehand similarly selects distant parcels.
Yet when I export the layer to non-virtual form (shapefile, spatialite table), the resulting layer works fine.
I've tried every combination I can think of for project and layer CRS because it feels like that's the likeliest cause, but I'm not fluent enough in projection and transformation to get a handle on it. Also I've searched high and low on this forum and others, to no avail.


texts:
I am volunteering for the wonderful Cambodian urban NGO Sahmakum Teang Tnaut as a GIS and data analysis advisor for 10 months. My Khmer counterparts have found that ArcGIS does not support Khmer Unicode for producing Khmer script in maps, such as for labels and titles. Does anyone have any experience with getting Khmer Unicode to work? I, and my colleagues have it installed on my computer so can use it in other programs, just not ArcGIS it seems.

At the moment we have to use any of the Limon typefaces, which are much more cumbersome to write in as it takes up to three characters per Khmer character.


texts: Some Edges of my Cropped Image Does not Provide the Entire Tile Based on the Shapefile

texts:
I have a shapefile and a big tiff image. I am trying to crop the tiff image based on the associated tiles of the shapefile. My problem is that I do not get the entire tile (tiles that are located in the edges of the tiff image)of the shapefile in my cropped images. It is a situation where the edge of the tiff image ends in the middle of a tile. In my code, I get the half that is in the picture but not the half in the tile. I always want to get the full.
Can you help me get the full tile?
Here is the code:
import os
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from shapely.geometry import mapping
import geopandas as gpd
import rasterio as rio
from rasterio.plot import plotting_extent
from rasterio.mask import mask
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
from pathlib import Path

shape_file = '/.../.geojson'
crop_extent = gpd.read_file(shape_file)
lidar_chm_path = '/.../reprojected.tif'
output = '/.../'

if not os.path.exists(output):
    os.mkdir(output)
i, j = 0, 0
while i < crop_extent.shape[0]:
    with rio.open(lidar_chm_path) as lidar_chm:
        try:
            lidar_chm_crop, lidar_chm_crop_meta = es.crop_image(lidar_chm,[crop_extent.iloc[i]['geometry']])
            lidar_chm_crop_affine = lidar_chm_crop_meta["transform"]
            j += 1
            lidar_chm_crop_meta.update({'transform': lidar_chm_crop_affine,
                       'height': lidar_chm_crop.shape[1],
                       'width': lidar_chm_crop.shape[2],
                        'nodata': -999.99})
            name = output + str(crop_extent.iloc[i][2]) + '_' + Path(lidar_chm_path).stem +'_' + str(j) + '.tif'
            with rio.open(name, 'w', **lidar_chm_crop_meta) as ff:
                ff.write(lidar_chm_crop[0], 1)
        except:
            None
    i += 1


texts: Handling multiple-extent-problem to create raster stack in R?

texts:
I am storing ~ 30 raster files in a folder and want to create a raster stack. I know that not all of these rasters have completely the same extent. They do have the same resolution and CRS.

Thus I try to set the same extent for all of my rasters in this folder with the following loop in R.

#reference raster to get the extent
a <- raster("F:/SDM/DATA/env_tiff/e_veg_d10.tif")

# where the new rasters will be stored
    outpath <- "F:/SDM/DATA/temp/"
    dir.create(outpath)

    files <- list.files(path="F:/SDM/DATA/env_tiff", pattern=".tif$")

    # add output directory
    outfiles <- paste0(outpath, files)
    # change extensions
    extension(outfiles) <- 'tif'

    for(i in 1:length(files)) {
      e <- extent(a)
      r <-raster(files[i])
      rc <- crop(r, e)
      rc<- mask(rc, e)
      rw <- writeRaster(rc, outfiles[i], overwrite=TRUE)
      print(outfiles[i])
    }

    #creates a raster stack with all environmental rasters in this file
    env_data <- list.files(path="F:/SDM/DATA/temp", pattern = 'tif$')
    env_data<- stack(env_data)

For my loop I get the following Error message:

  Error in (function (classes, fdef, mtable)  :    unable to find an
  inherited method for function ‘mask’ for signature ‘"RasterLayer",
  "Extent"’

When I comment the mask part, the loop works, but the extents are still different.

I am also not sure if I choose the right extent, since I don't know the minimum extent. I spent a lot of time searching a way to handle this multiple-extent-problem within my data and found the crop() and mask() functions.


texts:
I am trying to georeference a map from a Videogame to display it on a website, probably with openlayers.

The map in question is this Subnautica map.

I found a starting point here GeoServer create map from image. I managed to Reference the map in Qgis with some coordinates and export the layer to GeoServer and display it there.
My Problem is that I don't know how to proceed from there. I'd like to map the map 1:1 to "real world" coordinates, so that the center of this map is at 0,0 on the website. I think I have to chose a Map projection, like Web Mercator, and calculate where the coordinates from the game would be in the real world?

The map in the game is 4kmx4km in size.


texts: Arcgis 10 vs Arcgis 9.3 data compatibility

texts:
I have been doing some calculations such as slope in ArcGIS 9.3 but I needed to use ArcGIS 10 for fuzzy tool. My problem is when I open my data in ArcGIS 10, it shows me some wrong and extreme numbers like -3*10^38 to 3*10^38 while it is supposed to be 0 to 84.

What might be the problem?


texts: How to display the correct unit of measure in PyQgis?

texts:
Within a plugin I am developing I select features and for each line feature selected I get the feature distance using:

distance.measureLength(feature.geometry())

I also get the total length at the end using:

totallength = totallength + distance.measureLength(feature.geometry())

I display the total length in a message box but I don't know how to get the correct units of measure using pyqgis.  How do I know if I should use feet or meters?

I tried using:

distance.lengthUnits()

but that returns a value of 6.  Why?


texts: Create centerline from multi-digitized roadway lines Qgis 2.18

texts:
I am trying to create an absolute center line layer from US interstate data that is multi-digitized (center line for each direction).  I've looked at various options but they require a lot of processing involving voronoi polygons or adding points to lines and using a buffer or averaging those points.  They work on very small data sets but anything larger slows down considerably to the point where it would be faster to do it by hand.

I'm currently trying to get v.centerline in GRASS to work but the documentation on it is very limited so I'm not sure if it will accomplish what I'm doing.

It seems that Esri has a tool that will do this. Is there an equivalent in QGis?

The roads don't maintain a perfect parallel so any simple solution tends to fail in those areas.  Here's an example.


texts: Wrong results from TIGER geocoder for WA

texts:
For WA I'm getting wrong results from my local setup.
For the this address:

"1150 Varnum St NE, Washington, DC 20017"

SELECT g.rating, ST_X(geomout) As lon, ST_Y(geomout) As lat, (addy).* FROM geocode('1150 Varnum St NE, Washington, DC 20017', 1) As g;

I'm not getting any result:

SELECT g.rating, ST_X(geomout) As lon, ST_Y(geomout) As lat, (addy).* FROM geocode('1150 Varnum St NE, Washington, DC 20017', 1) As g;
 rating | lon | lat | address | predirabbrev | streetname | streettypeabbrev | postdirabbrev | internal | location | stateabbrev | zip | parsed
--------+-----+-----+---------+--------------+------------+------------------+---------------+----------+----------+-------------+-----+--------
(0 rows)

If I replace Washington, DC  with WA

1150 Varnum St NE, WA 20017

I do get a result:

SELECT g.rating, ST_X(geomout) As lon, ST_Y(geomout) As lat, (addy).* FROM geocode('1150 Varnum St NE, WA 20017', 1) As g;

rating |        lon        |       lat        | address | predirabbrev | streetname | streettypeabbrev | postdirabbrev | internal | location | stateabbrev |  zip  | parsed
--------+-------------------+------------------+---------+--------------+------------+------------------+---------------+----------+----------+-------------+-------+--------
     55 | -117.640723859113 | 48.0459212150471 |         |              | Vernon     | St               |               |          | Burien   | WA          | 98148 | t

But this result is actually wrong:

(-117.640723859113, 48.0459212150471)

Which are way different from the ones from the census site:

(-76.993, 38.94287)

I search if WA tables are loaded and got this result:

SELECT count(*) FROM tiger_data.wa_edges;
  count
---------
 2746400


texts: Using Check Geometries plugin for polygons in QGIS?

texts:
I've a few hundred polygon features with many duplicate nodes, gaps, overlaps and other geometry errors. Fixing all errors seems to be quite tricky. I had to run the duplicate nodes check twice to fix all duplicate nodes. After fixing minimum angles between segments (< 1,000000 deg) new duplicate nodes appear. The gap check isn't working because of a topology exception

  Input geom 0 is invalid: Self-intersection at or near point ...

but the self intersection check doesn't report any results. Fixing sliver polygons fails for most objects:

  Failed to merge with neighbour

What's the best way to use the 'Check Geometries' plugin?

It seems that tiny polygons may be highly problematic in some ways (e.g. polygons with three nodes containing duplicated nodes).

Should I start with a minimal polygon area check?


texts: Groupby Ploting - give each plot title/name

texts:
I have used this code to group by the continents and then plotted using geopandas. I would like to set the plot title for each continent. I would have liked also to draw the continents side by side.

import geopandas as gpd
import matplotlib.pyplot as plt
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

continent = world.groupby('continent')
continent.plot()
plt.show()


texts: Windows 10 GRASS init error on QGIS 3.10.4 install

texts:
I am a noob who is trying to reinstall QGIS 3.10 so that I can use some of the GRASS components but I keep getting this message:

GRASS init error : GRASS was not found in 'F:/QGIS 3.10/apps/grass/grass78' (GISBASE), provider and plugin will not work.

I am running Windows 10 with QGIS 3.10.4 with GRASS 7.8.2 and have reinstalled several times, I have gone into regedit and removed everything called qgis. I am at a loss but really could use the r.lidar components.


texts:
I'm using the GDAL SWIG python bindings to create a new GeoTiff raster for use in a blockwise exponential-decay function.

kernel_size = 177764054
driver = gdal.GetDriverByName('GTiff')
kernel_dataset = driver.Create(
    kernel_uri.encode('utf-8'), kernel_size, kernel_size, 1,
    gdal.GDT_Float32, options=[
        'TILED=YES', 'BLOCKXSIZE=1024',
        'BLOCKYSIZE=1024'])

# Make some kind of geotransform, it doesn't matter what but
# will make GIS libraries behave better if it's all defined
kernel_dataset.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS('NAD27')
kernel_dataset.SetProjection(srs.ExportToWkt())

kernel_band = kernel_dataset.GetRasterBand(1)
kernel_band.SetNoDataValue(-9999)

ERROR 1: TIFFNumberOfTiles:Integer overflow in TIFFNumberOfTiles
ERROR 1: GTiffDataset::Crystalize:No space for tile arrays
ERROR 1: TIFFNumberOfTiles:Integer overflow in TIFFNumberOfTiles
ERROR 1: TIFFReadDirectory:Cannot handle zero number of tiles
ERROR 1: TIFFNumberOfTiles:Integer overflow in TIFFNumberOfTiles
ERROR 1: TIFFWriteEncodedTile:No space for tile arrays
ERROR 1: /Users/jdouglass/Documents/habitat_quality_workspace/tmp/tmpvUVuqu, band 1: An error occured while writing a dirty block
ERROR 1: TIFFNumberOfTiles:Integer overflow in TIFFNumberOfTiles
ERROR 1: TIFFWriteEncodedTile:No space for tile arrays
...

Until finally:

More than 1000 errors or warnings have been reported. No more will be reported from now.

Admittedly, this is a very, very large raster (177,764,054^2 = 3.1600059e+16 pixels), which would also mean that there are a massive number of blocks (about 30,136,265,604).

Is the upper-bound of the number of blocks really 2^32 (4,294,967,296), represented by a 32-bit integer?  Is there any way around this apparent limitation other than to just use very, very large blocks?  Am I misusing my dataset creation options?

EDIT:

This raster is a convenient way to store the 2D kernel to a memory-efficient raster convolution (an input to pygeoprocessing.convolve_2d_uri).
The particular projection doesn't matter to me, but it is useful in that QGIS doesn't complain when the kernel (usually a kernel smaller than this!) is visualized so long as a projection is defined.
According to the GDAL GeoTiff format documentation, BigTiff=IF_NEEDED is supposed to be a default option in GDAL>=1.6.0, yet even when I specify BIGTIFF=YES as a creation option, I still get the same GDAL errors.


texts: Mapserver INCLUDE using multiple .map files

texts:
I am currently testing out a new way to display my Mapserver WMS's. Here is the context. I have a script that generates an individual .map file for each of my rasters. These mapfiles are fully functioning with their own Link, extent, layers, legends, etc... I'd like to combine two of the .map files into a "combined" .map file. I've been playing around with the INCLUDE call where I reference 2 of my generated .map files.

something like this:
MAP
..
INCLUDE 'Raster1.map'
INCLUDE 'Raster2.map'
...
END

This seems to work however mapserver is ONLY displaying the data for the first INCLUDE, NOT the second. Is there any way where I can get it to read both INCLUDEs?


texts: How do I generate a distance table using 2 different point files and a line file?

texts:
I have two point datasets, one containing multiple transit stops and the other containing multiple schools. I also have a street line dataset.

I need to generate a table that gives me the shortest distance from a school, to a transit stop, but only ALONG THE STREET LINE dataset. I know how to generate a "near table" but that gives me a straight line distance. That's all fine and good, but if I am a student who needs to walk to a transit stop, or drive to a large transit facility, a straight line distance does me no good. I need to be able to walk/drive along a street to reach my destination.

Keep in mind that I need to incorporate two different point files. A lot of the solutions only include one point file, just finding the distance between two of the points in the single data set. I need to use 2. The actual route for which the distance was calculated would be great too.

I am using ArcGIS 10.2 and have the Network Analyst Extension.


texts:
I have nine hillshades generated from DEMs and I was wondering if there is a way to merge all of them into a single layer. I tried using the Merge function but it didn't allow me select all the nine layers at once.

Is there any possibility of achieving this?

I am using ArcGIS Desktop 10.3.


texts: I would like to be able to create parallel lines with a defined distance, in the same way that the editor's parallel copy tool does.


texts:
I have an elevation raster (DEM) in PostGIS and want to compute a histogram using ST_Histogram. I use this query:

 SELECT (ST_Histogram( ST_Union(ST_Clip(rast,polygon.geom)),
         1, -- raster band number
         50, -- number of category/bins
         ARRAY[100] -- width category
        )).* as hist
 FROM elevation.dem,
      (SELECT ST_Transform(ST_SetSRID(ST_GeomfromText('Polygon((11 47,11  48 ,12  48 ,12  47 , 11 47))'), 4326 ), 3035) AS geom ) AS polygon
 WHERE ST_Intersects(rast,polygon.geom)

Result:

 min  | max  | count | percent
------+------+-------+-------
  478 |  578 |   782 |   5.7
  578 |  678 |  2533 |  18.6
  678 |  778 |  1489 |  10.9
  778 |  878 |   708 |   5.2
  878 |  978 |   635 |   4.7
  978 | 1078 |   635 |   4.7

However, I would like to have "pretty" categories like this:

 min  | max  | count | percent
------+------+-------+-------
  400 |  500 |   782 |   5.7
  500 |  600 |  2533 |  18.6
  600 |  700 |  1489 |  10.9
  700 |  800 |   708 |   5.2
  800 |  900 |   635 |   4.7
  900 | 1000 |   635 |   4.7

How can I do that?

SOLUTION

I found a quite hackerish solution with inspiration from Andre Silva. For my use case it seems to work:

SELECT min,
       max,
       round((percent*100)::numeric,1) as percent
FROM
(SELECT (ST_Histogram( rast,
         1, -- raster band number
         50, -- number of category/bins
         ARRAY[first_cat, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100] -- width category
        )).* as hist
 FROM
    (
    SELECT rast,
           ceil(((ST_SummaryStats(rast)).min)/100::float)*100 - (ST_SummaryStats(rast)).min AS first_cat
    FROM
        (
        SELECT ST_Union(ST_Clip(rast,polygon.geom)) AS rast
        FROM elevation.dem,
              (SELECT ST_Transform(ST_SetSRID(ST_GeomfromText('POLYGON((5.4931640625 45.07662388803723,8.32763671875 45.07662388803723,8.32763671875 42.76637205237248,5.4931640625 42.76637205237248,5.4931640625 45.07662388803723))'), 4326 ), 3035) AS geom ) AS polygon
        WHERE ST_Intersects(rast,polygon.geom)
        ) AS raster
    ) AS raster
) AS hg

Result:

 min  | max  | percent
------+------+---------
   -3 |    0 |     0.0
    0 |  100 |     3.4
  100 |  200 |     4.9
  200 |  300 |    10.5
  300 |  400 |     8.2
  400 |  500 |     6.7
  500 |  600 |     5.9
  600 |  700 |     5.0
  700 |  800 |     4.7
  800 |  900 |     4.8
  900 | 1000 |     4.6


texts:
In a virtual raster table's XML the LUT element can be used to specify a table lookup. To my understanding this could be used for some kind of reclass. This comes in really handy when different rasters are to be harmonized.

the Example Element in the  Section
0:0,2345.12:64,56789.5:128,2364753.02:255

However, I experience problems when the LUT section gets too long. In my case I want to put approx 1000 values in the list.
Does anybod know if there is a limit? And if so, which one?

If the limit is small, can the same raster be referenced more than once in a VRT and differen LUTs be appllied?


texts: Appending Data Into Existing Geometric Network

texts: I"m currently interning at a telecommunication company and have been tasked with converting their CAD maps into their GIS telecom model. I've been using FME to translate the CAD data with great success, but i've been roadblocked at the final step, which is appending the translated feature classes into the telecom model, necessary for the data to be uploaded to our Oracle database. The data that I'm outputting from FME is in simple line feature format and i'm trying to append these features to their corresponding empty feature classes that are in complex edge format. If you have any ideas please let me know.


texts:
I have two grids.

One, call it A, is 4320x2160 cells and spans the whole globe.

The other, call it B is 8640x3432 and spans the circumference of the globe, but is clipped toward mid-latitudes.

Both files are in ArcGrid ASCII format. The clipped file has a header of:

ncols         8640
nrows         3432
xllcorner     -180
yllcorner     -58
cellsize      0.0416666666667

Looking at the grids in QGis, I see that the B contains four cells for every cell of A.

I would like to sum each group of four cells in B to create a new grid that overlays directly with A and has a one-to-one correspondence of cells.

Is there a way to do this in Linux using either Python or command line only tools?


texts: Putting scale bar labels below scale bar in QGIS?

texts:
I added two scale elements to my map composition. For layout harmonisation, I want the ticks and the label text below the scale bar.

Standard is this:

What I want is this (manually changed in the image):

Is there an option in QGIS for this?


texts: How to calculate area of intersection between raster and polygon in a web map?

texts:
I need to calculate the area of intersection between a vegetation raster and a polygon drawn by the user in an ArcGIS Server JavaScript API web-map, as in this mockup:

After the user has drawn the polygon, the web map should return the total area of vegetation within the polygon.

I'm in the early stages of investigating this problem so I don't yet know the final data formats, but assume that the raster contains 1 to show where vegetation exists and NODATA where it doesn't, and that is has been projected into square metres.

Presumably the raster will be served as a tiled layer, although I'm open to other suggestions (Image Server?).

I don't believe there is any functionality in the ArcGIS Server JS API to calculate the area of intersection between a raster and a vector. In that case, is the correct approach to write a geoprocessing function to handle this on the server?


texts: How to split a (road) buffer using the centerline?

texts:
I have a road centreline dataset which I have buffered using the road width. I used the standard QGIS buffer tool. I realise for my analysis I actually need to split the road buffer in half and essentially create a 'right' and 'left' side of the road.

My idea is to use the QGIS processing toolbox buffer and select the 'don't make caps' option to create a square buffer. Then I should be able to use the original centreline to split the buffer.

My questions are, if anyone has done something like this:

A. Is there a plugin that can already do this? Or if not,
B. What would be the best way to run this process over 30,000 polylines? And,
C. How could I label each side of the buffer? (keeping in mind the labeling could be as arbitrary and 'side A' and 'side B'.)

I'm also open to using ArcGIS or any other tool if it means I could get this done faster.


texts: Create New Feature Class Based off Multiple Feature Classes

texts:
I've searched previous posts to see if anyone asks a similar question, but I haven't found anything. I have a multiple feature classes, each containing multiple polygons. I have two desired outcomes:

A) I would like to find any easy way to create a new single feature class in areas where more than 50% of my input layers intersect. The number of input layers I have will vary, but let's use 6 as an example. It would be nice to be able to create a single feature class no matter which of my four input layers make up the more than 50%. I've tried to create a model using Intersect, but I am missing the logic piece that tells the system that not all 6 inputs need to overlap, just at least four.

B) The other output I would like to have is a new feature class which creates new polygons which are an average area based off my inputs. In this case too I would like this tool to create a new feature class with the mathematical average of areas where more than 50% of my inputs overlap. I've searched, but cannot find any such tool that can get me at my output.

I currently have an ArcInfo (Advanced) license, but no extensions. If there is a way to accomplish my goals with either ArcView (Basic desktop) or ArcInfo that would be preferred. But if there is an extension out there I would be curious to know what it is.


texts:
I would need some help with the Reclassify tool. I have written a script that produces rasters which I need to reclassify. My problem is, I want the values to be reclassified by standard deviation which I don't know how to do. I would like to make two new classes:

values < mean + 2 * standard deviation

values >= mean + 2 * standard deviation

I need help with putting this into a script.


texts: Is it possible to rotate items in QGIS Print Composer?

texts:
Is it possible to rotate items in Quantum GIS Print Composer?

I just didn't found and before submit feature demand just double-check this.


texts: OSM: getting land area polygon of a country/region

texts:
I was able to get the polygon of a country and its subdivisions using the Overpass API, overpy Python package, and ISO codes with the following queries:
(rel[boundary="administrative"]["ISO3166-1:alpha2"="DK"][admin_level="2"];>;);way(r);(._;>;);out;

(rel[boundary="administrative"]["ISO3166-2"~"^DK"];>;);way(r);(._;>;);out;

However, in some cases, the polygon includes the maritime areas as well. I would like to retrieve only the land area.
I'm aware of the land_area tag but cannot figure out how to use it along with the ISO3166-1:alpha2 tag without receiving empty results (is it even possible)?
Any solution to get the land area polygon of a particular country/region? I know I can download the polygons on some web sites but I need to implement a solution in a Python script preferably...
The situation seems to change according to the country (I guess the tags are not consistent?) E.g. In the UK, Norway, Sweden, the maritime area is included at both, country and regional levels (in brown below). In France, Ireland, Denmark, the maritime area is included at country level (in green below) but not at regional level.


texts: Query a different density of points based on zoom level in PostGIS

texts:
I have a database (using PostGIS) of a few million points that I load onto an interactive map. Obviously the map can't handle this many points, so I currently restrict it by limiting to the points closest to the center of the map. This works great when zoomed in, but shows a tiny cluster of points when zoomed out.

Instead, I'd like to show a sampling of points that represents the visible area.

Here's an example of the functionality I want to copy. As you zoom out, points that are nearby are merged into a single point and as you move in, you can see more detail.

https://www.lds.org/maps/meetinghouses/

How can I accomplish this in a performant way?


texts:
I am trying to get the default marker NOT to show, but it's coming up WITH the icon that I want, so I have a marker and icon like this:

I'm getting my data from a OSM nomination search result.

if (response) {
        var mArray = [];
        for (var i = 0; i < response.length; i++) {
            var result = response[i];
            var s = result.display_name;
            var ex = result.extratags;
            var lk = ex.website;
            var ph = ex.phone;
            var hr = ex.opening_hours;
            var poly = result.geojson;
            if (hr) {
                var hours = hr.replace(/\;/g, '
');
            }
            var arr = s.split(', ');
            var dd = s.substring(0, s.lastIndexOf(", "));
            var cc = dd.substring(dd.indexOf(",") + 1);
            var rIcon = result.icon;
            var myIcon;
            if (rIcon) {
                myIcon = new L.Icon({
                    iconSize: [35, 35],
                    iconAnchor: [13, 27],
                    popupAnchor: [1, -24],
                    iconUrl: rIcon
                });
            } else if (!(rIcon)) {
                myIcon = new L.Icon({
                    iconSize: [35, 35],
                    iconAnchor: [13, 27],
                    popupAnchor: [1, -24],
                    iconUrl: 'lib/images/question.png'
                });
            }
            var srchPops = "
";
            srchPops += "

" + arr[0] + "
";
            srchPops += "

" + cc +
                    "

";
            srchPops += "<button class='btn btn-default btn-xs add-field' value='" + result.lat + ", " + result.lon + "'>Add Point</button>";
            srchMarkers = new L.marker([result.lat, result.lon], {icon: myIcon}).bindPopup(srchPops);
            srchPoly = L.geoJson(poly);
            mArray.push(srchMarkers, srchPoly);
        }
    }
    srchresultMarkers = L.layerGroup(mArray);

So what am I doing wrong?


texts: Updating R version in R-ArcGIS/r-bridge?

texts:
I recently installed the R-ArcGIS bridge, following the steps on the associated github. I have ArcGIS 10.3, along with both R 3.2 & R 3.3 installed. It is only recognizing the R v3.2 installation, but I need it to recognize R v3.3. I have tried using the 'Update R bindings' Arc tool that comes in the installation, but it only recognizes version 3.2. The current R-version in the registry is 3.3.

Is R-bridge compatible with the newest versions of R?

Has anyone experienced issues changing the associated R version with R-bridge?


texts: Determining Min and Max elevation in Python using GDAL?

texts:
I have a DEM file and shapefile of specific area, and I want to
Create a module in Python that writes name of each area which is shown in the picture, with the greatest differences in altitude on the territory (min-max range of the altitude).

I also have to create a new layer containing the shape of these areas, together with all its original attributes.

I have GDAL and ArcPy libraries.

Could you please help me determine what I should do?


texts:
I'm attempting to complete a batch mosaic operation in ArcMap 10.4.1 Python window. I need to run the Mosaic tool about 80 times for sets of three rasters. Basically I want the second and third raster to mosaic to the first raster with all of the default mosaic properties. My code to do this is pasted below, but after playing around with the inputs for the Mosaic Management tool I keep getting back a parsing error SyntaxError in the line where the Mosaic operation occurs.

# Name: Batch Extract by Mask

... # Description: Iteratively calculate extracts the flow velocity of cells
... # Requires spatial analyst
... # MAKE SURE RASTERS FOR Connectivity ARE LOADED IN THE MXD TABLE OF CONTENTS
...
... # Import system modules
... import arcpy
... from arcpy import env
... from arcpy.sa import *
... import math

... # Set environment settings

... env.workspace = "C:\Subbasin Connectivity Flow Time"

... # Check the ArcGIS Spatial Analyst extension license

... arcpy.CheckOutExtension("Spatial")

... count = 1

... # Run for each raster in the directory

... while count < 80:
... # Execute times flow direction calculation
...     name1 = "Day72Subbasin" + str(count) + " Bin 1.tif"
...     Bin1 = Raster(name1)
...     name2 = "Day72Subbasin" + str(count) + " Bin 2.tif"
...     Bin2 = Raster(name2)
...     name3 = "Day72Subbasin" + str(count) + " Bin 3.tif"
...     Bin3 = Raster(name3)
...     arcpy.Mosaic_management(Bin2;Bin3,Bin1,"LAST","FIRST","0", "", "NONE", "", "")
...     count = count + 1


texts: Iterating over layer features, zooming fully on map and exporting them as PNG images in QGIS

texts:
I have 2 layers in QGIS map, One is raster layer and other is a vector layer having location of sites. Both are having same CRS. I want to zoom in to each feature in the vector layer and save the map. I am using QGIS 3.8
I also tried the QTimer function as well as indicated on How to loop through features in a layer, zoom to that feature and export to png
but it doesn't seem to be working.

Here is the code

    ids =None
    def get_screenshots(layer):
        global id
        ids = layer.allFeatureIds()
        def savemap():
            qgis.utils.iface.mapCanvas().refreshAllLayers()
            layer.select(id)
            qgis.utils.iface.mapCanvas().zoomScale(400)
            qgis.utils.iface.actionZoomToSelected().trigger()
            qgis.utils.iface.mapCanvas().saveAsImage("C:/Users/sumitsingh/
            Desktop/screenshots/pic_{}.png".format(id)
            layer.deselect(id)
        for id in ids:
            savemap()

    l = iface.activeLayer()
    get_screenshots(l)


texts: ArcGIS Pro Export Map PDF with extra colored lines

texts:
I am exporting map to pdf using ArcGIS Pro at 600 dpi resolution. The exported map is coming with some thin colored lines.
I tried to export same map to png it appeared fine.

Am i missing something here regarding export settings for pdf ?

I have an interpolated raster these lines are reflecting the colors of that raster.

When i un-check that raster then it prints pdf fine.

I cannot include complete map screen shot but i am including the area where i dont have data. See the green lines in image.

When i select these blocks and right click i get the option to copy image. Seems like different small images have been stitched to generate this pdf.


texts:
Performing research on a new Python script to use in QGIS where I need to create a random point based on another point. This new point will need to have its location based on a random degree (two user inputs with a range from 0 to 360 degrees), as well as a random distance in meters from the original point (two user inputs from 1 to 100 meters).  The result should be a new layer or appending an existing layer.

Been finding some snippets of code I think I can start with but searched pyQGIS cookbook and python in general but so far no luck regarding random degrees.

Is there some sort of function for a random degrees selection?


texts: Some of the CartoDB maps online, like Illustreets react constantly to mouse hover instead of requiring a click. Is there a way that I can do that with my CartoDB maps?


texts: Defining partial or complete overlap with QgsSpatialIndex

texts:
I'm working on adding some functionality to a QGIS plugin. It currently uses a QgsSpatialIndex to efficiently compute whether the features of a vector layer intersect (potentially millions) of QgsRectangles representing tile bounds. This is reasonably straightforward:

spatialIndex = QgsSpatialIndex()
canvas = self.iface.mapCanvas()
geomTransform = QgsCoordinateTransform(canvas.mapRenderer().destinationCrs(),\
                QgsCoordinateReferenceSystem('EPSG:4326'))
# layer is a QgsVectorLayer, from a combobox selection
vals = {feature.id(): feature for (Feature) in layer.getFeatures()}.values()
for v in vals:
    v.geometry().transform(geomTransform)
map(spatialIndex.insertFeature, vals)

Then later I can check whether a given QgsRectangle intersects features participating in this index:

# includes pseudo-code
if spatialIndex.intersects(rectangle):
    # Feature intersects the rectangle
    # Both partial and complete overlaps are included
    if complete_overlap:
        render_one_way()
    elif partial_overlap:
        render_another_way()
else:
    # No features intersect the rectangle
    dont_render()

This is very fast. However, I need to be able to distinguish between QgsRectangle objects that intersect features in the spatial index partially, and those that overlap completely. The QgsSpatialIndex.intersects method only tells me if there is any kind of intersection.

How can I efficiently determine whether it is a partial or complete overlap?


texts:
I'm trying to install a version of the QGIS API that is compatible with the python application I'm developing in Anaconda 2 (Python version 2.7). When I install qgis using conda, it indicates that the QGIS API is only compatible with Python 3. I realize that I eventually need to convert my code to Python 3, but I'm looking for a short term fix. Is it possible to install a version of the QGIS API via conda that is compatible with Python 2.7? When I use the command conda install -c conda-forge qgis=2.18.22 or similar it doesn't find anything.

thanks.
Jack Andrewson


texts: Unresponsive Contents pane in ArcGIS Pro

texts:
I’m having an issue where the new install of ArcGIS Pro version 2.7 isn’t responsive in certain areas. I’m using Windows 10, and my machine could handle the workload before upgrading.
Basically the Contents pane will continuously give me the loading wheel and be unable to interact with. Most buttons on the top ribbons are grayed out as well. The only things I can do are panning the map and geoprocessing tools - which are considerably faster than before. Which is good for some situations, and makes others impossible. The other GIS analyst who sits next to me also has this exact issue.