texts: GRASS within QGIS "Python error: An error has occurred"?

texts:
I'm very new to using GRASS within QGIS. I want to buffer a vector layer. I want to use the GRASS GUI as I can use the -c option to not create caps however I cannot launch v.buffer.distance.

I'm guessing I need to re-install something? I've just done a clean install of QGIS 2.16.0


texts: Get GeoJSON data from AJAX enabled service

texts:
I want to get GeoJSON data from GeoServer using asp.net and AJAX enabled service. firstly I made a class for request web:

public class Request
{
    public static string WebRequestByUrl(string url)
    {
        string jsonData = "";
        try
        {
            var request = WebRequest.Create(url);
            request.ContentType = "application/json; charset=utf-8";
            var response = (HttpWebResponse)request.GetResponse();
            using (var sr = new StreamReader(response.GetResponseStream()))
            {
                jsonData = sr.ReadToEnd().Replace("'", "-");
            }
        }
        catch (Exception)
        { }
        return jsonData;
    }
}

and AJAX enabled service codes

[ServiceContract(Namespace = "")]
[AspNetCompatibilityRequirements(RequirementsMode = AspNetCompatibilityRequirementsMode.Allowed)]
public class Mapping
{
    // To use HTTP GET, add [WebGet] attribute. (Default ResponseFormat is WebMessageFormat.Json)
    // To create an operation that returns XML,
    //     add [WebGet(ResponseFormat=WebMessageFormat.Xml)],
    //     and include the following line in the operation body:
    //         WebOperationContext.Current.OutgoingResponse.ContentType = "text/xml";
    [OperationContract]
    public string WebRequestByUrl(string url)
    {
        // Add your operation implementation here
        return Request.WebRequestByUrl(url);
    }

    // Add more operations here and mark them with [OperationContract]
}

and JavaScript side

var mapping = {
getTmJson : function () {
    var url = "http://10.212.81.15:8080/geoserver/BEIS_WS/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=BEIS_WS:ilce&maxFeatures=50&outputFormat=application/json";
    var service = new Mapping();
    service.WebRequestByUrl(url, function (data) {
        console.log(data)
        mapping.showData(data);
    }, null, null);
},
showData: function (data) {
    var raster = new ol.layer.Tile({
        source: new ol.source.MapQuest({
            layer: 'sat'
        })
    });

    var vector = new ol.layer.Vector({
        source: new ol.source.GeoJSON({
            projection: 'EPSG:5254',
            object:data
        })
    });

    var select = new ol.interaction.Select();

    var modify = new ol.interaction.Modify({
        features: select.getFeatures()
    });

    var map = new ol.Map({
        interactions: ol.interaction.defaults().extend([select, modify]),
        layers: [new ol.layer.Tile({ source: new ol.source.MapQuest({ layer: 'osm' }) }), vector],
        target: 'map',
        view: new ol.View({
            center: [0, 0],
            zoom: 2
        })
    });
}

}

and result I can get JSON data from GeoServer, but I can't make it integrate with OpenLayers vector

var vector = new ol.layer.Vector({
        source: new ol.source.GeoJSON({
            projection: 'EPSG:5254',
            **object:data**
        })
    });

error is
** Uncaught TypeError: Cannot read property 'c' of undefined**

where is my mistake?


texts: QGIS - Combining contours without elevation and point elevation data

texts:
I have a geopackage (Ordnance Survey Terrain 50) that has two layers: a layer with contours but no elevation data, and a layer of point elevation data.

I'd like to create a DEM. The technique I've found requires contours that have elevation data. How do I combine the two layers in QGIS to give me this?


texts: GeoServer - set declared SRS with rest api

texts: I am using a c# application that creates and update featuretypes in Geoserver via REST. I need to set the declared srs, but I couldn't find how to do so.


texts: Change attribute values of multiple fields via PyQGIS

texts:
I am starting with PyQGIS and I would like to change the values of two fields within selected features:

from qgis.PyQt.QtCore import QVariant
parent = iface.mainWindow()
mc = iface.mapCanvas()
myLyr = mc.currentLayer() #aktueller Layer, der als Quelle genutzt wird
flds =myLyr.fields()
fld1 = "BTYP" #Name of the first field
fld1i = myLyr.fields().indexOf(fld1) #get the index of the first field
fld1val = "herdamit"#Value of the first field
fld2 = "WS"
fld2i = myLyr.fields().indexOf(fld2)
fld2val = "superduber"
caps=myLyr.dataProvider().capabilities()
myData = {fld1i: fld1val,fld2i: fld2val}
mySel = myLyr.selectedFeatures()
print(myData)
with edit(myLyr):
    if caps & QgsVectorDataProvider.ChangeAttributeValues:
        for f in mySel:
            myLyr.dataProvider().changeAttributeValues({f : myData})
            myLyr.updateFeature(f)


texts:
I am hoping for an answer here even though it's a long shot. My goal is to add a raster service into ArcMap that shows the Doppler RADAR information that is shown at the Weather Office website.

I am hoping that this data can be dynamic, so that I can do analysis of it and so that it refreshes every ten minutes or so.

I work for a Government organization, so there should be no private sector issues.

Over the last month I've sent out quite a few emails to different possible data sources. Each one has not really had an idea about what I am inquiring about.

Here is an image of what I am looking for. If you want to check out the website click here.


texts: More efficient approach for processing large polygon dataset

texts:
Using ArcGIS Desktop 10.6, I have the following workflow in a python script to try and calculate coverage area affected on a wireless network if certain towers went down:

User enters a list of polygon IDs (these are the tower outages)
Using an in-memory layer, script selects all polygons from a layer in a
file geodatabase, except the ones listed by the user

Important note: Polygon layer is 500GB in size and cannot be generalized.
 Contains complex geometries and often intentional overlap between polygons

Polygons from #2 are merged (arcpy.Merge_management)
Loop through all the polygons that were removed and clip any overlapping areas from merged result in #3
Merge output from #4 and calculate its area.

I am fairly confident that I have optimized the script as much as possible, and due to the size of the data, it is taking a significantly large time to process a result.

The question I have is what is a good alternative approach that will significantly reduce the processing time for the user but still accomplish the same workflow?

Things I am currently contemplating:

Pre-cooking every combination of result using a union on the original layer.  This may be a beast to process, but would only need to happen every ~6mths.
Storing in SQL Geometry and using native SQL spatial functions to return results.  Not tested, but I assume some of the more complex geoms might cause some issues.
Raster based approach?  Not put a lot of thought into this one.


texts: Is it possible to create a read-only project in QGIS? Or is there any viewer software available for QGIS project file? What I am after is a read-only project that shows and keep labels and rendering.


texts: The chart takes much time to update with the ClickMap function in GEE, is there any way to optimize this process?

texts:
I'm working in a code to analyze NO2 emissions in Mexico. I had used the clickMap function in GEE to update the chart with data of the state I choose. My issue is that it takes so much time to update the chart. I don't know if there is something wrong in my code, here it is: https://code.earthengine.google.com/da23c5d31ca46d705cb25f1c47066311.
So far, I tried to change the scale but still it's slow. Is there any way to make it faster?

Thanks in advance for your comments!


texts:
I am new to ArcPy and was trying to get a new layer from two shapefiles. The new layer would only have the data if the two layers intersect and one of the attributes ("VEGETATION") = "NV".
import arcpy
import numpy

#Locations
arcpy.env.workspace = "c:\Users\C\Desktop\GIS\Geodatabases"
LClocation = "c:\Users\C\Desktop\GIS\Geodatabases\ex_LandCover.shp"
RBveglocation = "c:\Users\C\Desktop\GIS\Geodatabases\ex_RedBluffVegetation.shp"

a_RBVeg = arcpy.MakeFeatureLayer_management(RBveglocation, "VEGETATION")
a_LC = arcpy.MakeFeatureLayer_management(LClocation, "LU_CODE")
updateveg = arcpy.SelectLayerByLocation_management(a_RBVeg,'INTERSECT', a_LC, 0, 'NEW_SELECTION')

query = '"VEGETETATION" = "NV"'
updatevegQ = arcpy.SelectLayerByAttribute_management(a_RBVeg, "SUBSET_SELECTION", query)

arcpy.CopyFeatures_management(updatevegQ,'updateveg4')


texts: How to save a GeoJSON layer with the select option for the items in a field in QGIS

texts: I have a GeoJSON layers with many items. If you open the attribute table you can select an option in one field clicking over the field and the choosing one of the displayed options. The problem is that  I don't know how to save the GeoJSON  layer with the options that I chose. I want to save the options chosen so when I open the layer again, the chosen options will be the selected ones by default.


texts: Reclassifying raster stack based on condition and other layers using R?

texts:
I want to reclassify raster r1 in the following form:

all the values that are larger than r2 to be 1
all the values that are less than r3 to be 0
all the rest of the values to be equal to r4

I use overlay to set the values for the first two conditions, but I can't do that for the third one. I am also interested to know how to put all these in one function.

library(raster)
r1 <- raster(nrow=5, ncol=5)
r1 <- setValues(r1, runif(ncell(r1)))
r2 <- setValues(r1, runif(25 ,0.6,0.9))
r3 <- setValues(r1, runif(25 ,0.2,0.4))
r4 <- setValues(r1, runif(25 ,1,2))

x <- overlay(r1, r2, fun=function(x,y){ x[x>y] <- 1; x})
x2 <- overlay(x, r3, fun=function(x,y){ x[x<y] <- 0; x})


texts: What equivalent zoom level would ISS astronauts would be seeing the Earth?

texts: Like for Google Maps, there are 20 zoom levels. What zoom level would approximately give us the view that astronauts have?


texts: Is there a command line tool like `gdalinfo` for vector files?

texts: I use gdalinfo all the time to obtain basic information about raster files. I was wondering if there was a similar tool to obtain information from vector files, such as the spatial projection, list of attributes and units, the types of geometries, etc? I can of course put a vector file into QGIS or use geopandas, but was just wondering if some command line tool already exists that is the equivalent of gdalinfo. I suppose I would prefer something in python if available.


texts: How to calculate distance from lat/long point to each habitat category in QGIS?

texts:
I am very new to QGIS and am struggling so hopefully someone will be able to help me with what I assume is a relatively simple task.

I have a number of vector layers of land cover from the land cover 2007 dataset, available from edina digimap. (categories: agriculture, alternative, builtup, freshwater, grassland, intertidal, offshore marine). I also have a vector layer of points (lat/long). I basically need to calculate the distances (km) from each point to each habitat category, and end up with a table with column names (lat, long, dist to hab A, dist to hab B, dist to hab C....etc.).

I've been trying to convert the landcover shapefiles to rasters as this seems to be the first step (from looking online), and have thus far been unsuccessful - I think because I'm not sure what values to use to define extent and the like. Having trawled online (similar) questions I then tried to convert to rasters in R (more familiar with this environment) and couldnt load the files in - fell at the first hurdle!

Hoping someone can help me out!


texts: When I try to use the nvis module in Grass, QGis, pops-up a warning window saying "Module nviz not found". how can I solve this issue and use the nviz?


texts: Check geometries in QGIS for multiple KML Files all at once in R or Python?

texts:
I am very new to GIS analysis so please bear with me.

I have 2000 KML files.

I am currently using QGIS Check Geometries function to generate the error report that consists of the following.

A CSV File
A Geopackage file
A Shapefile

Is it possible to do this in R or in Python for all 2000 KML files at once?


texts: Automating select polygons by distance within same layer using QGIS?

texts:
I have say, 6,000 polygons roughly 25 acres each, all touching one another, basically a grid but with irregular shapes.

Using QGIS, how do I start with one shape then in an automated process, select only those shapes that aren't within 300 feet?

So of all 6,000 shapes how many can be "used" based on the fact they aren't within 300 feet of each other.


texts: Making Images in Collection Binary Based on Threshold

texts:
I am trying to make every image in a collection into a binary image based on whether it's band's pixel value is higher or lower than a certain threshold. For instance, if an image's pixel in a selected band is above 450 I want it's value to be reassigned to 1, but if it is below 450 I want it to equal zero. The code I wrote is able to make all pixels below 450 equal to 0, but I can't get it to do the conversion to a pixel value of 1 for the greater than portion...
Here is the code I wrote:
https://code.earthengine.google.com/1cc75cf12cd5891bb0da252344018f01


texts: Distance to nearest hub output not correct

texts:
In QGIS 2.6, I'm working with two point layers: center of the population per census tract and stores that accept food stamps. I'm using "distance to nearest hub" to calculate nearest food stamp outlet from the population center of each tract. I'm using miles as my measurement. However, looking at the output attribute table after running the analysis, distances appear to be very small- something like 0.00003. This obviously isn't right. Using the measure tool, I measured distances between two points in miles, then I went to google maps to measure the distance between the same two points, and the miles are close.

I thought it could be with my CRS, but everything lines up correctly, and I have CRS on the fly enabled. Specifically, I'm using WSG 84 (EPSG 4326) for all layers. Any ideas on this issue?

UPDATE: Pictures added based on responses below.


texts: Exporting shapefile/feature class polygon to XY pairs/ASCII

texts:
For a polygon shapefile, I'm needing a polygon vertex coordinate output file like this for a legacy program:

polygon1, AUTO
  -120.3434    34.694354
  -120.3750    34.694354
  -120.3750    34.713433
  -120.3434    34.713433
END
polygon2, AUTO
  -121.3434    35.694354
  -121.3750    35.694354
  -121.3750    35.713433
  -121.3434    35.713433
END

When I try Spatial Statistics Tool's Export Feature Attribute to ASCII, however, the output looks like output for a point file, ie just one XY pair for each polygon (each polygon has many vertices). In another program I could add in the , AUTO and END pieces - I'm needing most just how to export the XY polygon pairs.

I've come across an Ungenerate tool that might be another starting point, but ArcMap 10.6 won't seem to open it from the tool search display, and also a Feature Vertices to Points tool, but looks like that requires the Advanced license.

Any ideas?


texts: After using Apply Symbology From Layer in ArcGIS Pro layer becomes invisible

texts:
The layer still exists, i can open the attribute table, zoom to the layer etc.
And if i change the Symbology manually, it reappears.
Anyone knows why and how to fix it?


texts: Capturing Building level information (name and street names) from Basemap Using QGIS?

texts:
I wanted to capture House number and street number of Malaysia. I Found that Heremaps  streets and Open route services have proper House numbers and street names to capture like this in the below image

Can we take those services to QGIS as basemaps? Is there any WMS  or any services that can help to capture the information? Is there any plugin that can help you to download the information from WMS?

I have tried this GDAL for Heremaps  Heremaps as basemap, but i am unable to add the Heremaps Terrain to QGIS. It would be very helpful if i can add the heremaps Terrain as it is to the QGIS as a basemap.


texts:
I made a python plugin for qgis 3.4 and I built the interface in Qt design.
I have the problem when I insert an image.
I create a Label in qt design and configure the directory on pixmap, on Qt design it is displayed well.
I perform the refresh on qgis but when I open the plugin the image cannot be seen.
Who can help me?


texts: Corrupted machine, unable to find "Recover Lost license"

texts: My machine I use for Arc Desktop recently... crashed. I am trying to authorize my license on a new machine. This is a ArcGIS Desktop for Personal Use license (single use). Below is a screenshot of the licensing tab (with no option to recover lost licenses). Any help?


texts: Piping one SAGA Raster Calculator output into another chaining algorithms in QGIS processing script?

texts:
I am using QGIS 2.18.15.

I want to streamline multiple SAGA raster calculator calculations using the output of one calculation (not loading it into the table of contents) as input raster for the next and so one. I found some useful things at Pipe output of one processing algorithm as input into another algorithm in QGIS 2.18 Processing script and in the QGIS doc here.

I built a script looking like this:

##MaechtigkeitunterhalbeinerHoehe=name
##Input_Raster=raster
##Hoehe= string 28.65
##Pos_Differenzen=output raster

formula='a-{}'.format(Hoehe)

Differenzen=processing.runalg('saga:rastercalculator',Input_Raster,[],formula,0,False,7,None)

Pos_Differenzenalg=processing.runalg('saga:rastercalculator', Differenzen['Result'],[],'a*(a>0)',0,False,7,Pos_Differenzen)

As you can see I am trying to access the output of "Differenzen" by using Differenzen['Result'] as the input parameter for "Pos_Differenzen". When I define an output rater parameter like ## Differenzenout=output rasterit at least generates no error through the first part, but it load the result into the table of content (which I am trying to avoid) and also gives me an error for the second part. I have also tried 'OUTPUT' instead of 'RESULT'but to no avail. If I run the second algorithm on its own with the previously calculated temporary raster from the Diffrenzen algorithm as ##Input_raster it works fine again.

What am I doing wrong?


texts: software architecture for creating routes?

texts:
I have a system that can alert staff by email if their device goes into/out of a geo-zone

However, we now need to create a route, and then email if someone goes off this route

Technically, how would i do this? I have no idea!
Do you just check their GPS location and see if it exists within the GPS co-ordinates of a route?

I'm using openlayers framework

....an afterthought
Is a route a really narrow geozone?


texts: Showing QGIS map in Qt/C++ custom widget with QgsMapCanvas?

texts:
I'm trying to add a QGis map in Qt/C++ widget using QgsMapCanvas, but I have some problems...

I'm working with the following environment:
- Ubuntu 12.04
- QGis 2.8.1
- Qt Creator 2.4.1 (based on Qt 4.8.0 - 32 bit)

I found online a sample code that should do what I'm looking for, but something doesn't work..

The code is the following:

int main(int argc, char **argv)
{
    QgsApplication app(argc, argv, true);

    QString myPluginsDir        = "/home/pesce/qgis";
    QString myLayerPath         = "/home/pesce/qgis_sample_data/shapefiles/alaska.shp";
    QString myLayerBaseName     = "al";
    QString myProviderName      = "ogr";

    // Instantiate Provider Registry
    QgsProviderRegistry * preg = QgsProviderRegistry::instance(myPluginsDir);

    QgsVectorLayer *mypLayer =  new QgsVectorLayer(myLayerPath,myLayerBaseName,myProviderName);

    QgsSingleSymbolRendererV2 *mypRenderer = new QgsSingleSymbolRendererV2(QgsSymbolV2::defaultSymbol(mypLayer->geometryType()));

    QList <QgsMapCanvasLayer> myLayerSet;

    mypLayer->setRendererV2(mypRenderer);

    if(mypLayer->isValid()){
        qDebug("Layer is valid");
    } else {
        qDebug("Layer is NOT valid");
    }

    // Add the Vector Layer to the Layer Registry
    QgsMapLayerRegistry::instance()->addMapLayer(mypLayer, TRUE);
    // Add the Layer to the Layer Set
    myLayerSet.append(QgsMapCanvasLayer(mypLayer, TRUE));

    // Create the Map Canvas
    QgsMapCanvas * mypMapCanvas = new QgsMapCanvas(0, 0);
    mypMapCanvas->setExtent(mypLayer->extent());
    mypMapCanvas->enableAntiAliasing(true);
    mypMapCanvas->setCanvasColor(QColor(255, 255, 255));
    mypMapCanvas->freeze(false);
    // Set the Map Canvas Layer Set
    mypMapCanvas->setLayerSet(myLayerSet);
    mypMapCanvas->setVisible(true);
    mypMapCanvas->refresh();

    return app.exec();
}

The code is correctly compiled and the program correctly starts, but the map is empty..nothing is show... the selected layer (/home/pesce/qgis_sample_data/shapefiles/alaska.shp) seems to be not valid (mypLayer->isValid() returns "false")... but it's a .shp file....it should be valid!

Which kind of files are supported?


texts: How to mosaic/align multispectral georeferenced images into one merged file with GDAL

texts:
I am trying to work my head around how to make mosaic from georeferenced multispectral images. Header of those images already contains information about GPS location.
Separate image is saved as *.tif file and consist of 5 separate bands. Ultimately I would like to make one big TIFF file which would include all multispectral images merged into one. Also I would like to execute merging algorithmically with the use of Python.
For purpose of testing I am just trying to align two of those multispectral images (pic1.tif and pic2.tif) into one. See a screenshot of pics: pics.

With some googling I found out that this can be done with the use of GDAL. Since my images already include information about location I just excecuted
command = gdal_merge.py -o mosaic.tif -of gtiff pic1.tif pic2.tif
output = subprocess.check_output(command)

The command executed succesfully, but the output image (mosaic.tif) sems to be equal to pic2.tif.
After that I also tried to manualy insert GCP points with:
gdal_translate -gcp 620.5 456.0 33.98 44.839 -of GTiff pic1.tif newpic1.tif
gdal_translate -gcp  620.5 456.0 34 44.63 -of GTiff pic2.tif newpic2.tif

where the numbers: 620.5 - 456.0 - x - x represents respectively: image_center_x - image_center_y - longitude - latitude
After execution of:
python C:/pathtofile/gdal_merge.py -of GTiff -o mosaic.tif newpic1.tif newpic2.tif

also no luck.
Can tell me what I am doing wrong and what should I change/do, to make mosaic?
I am open to another solution if GDAL is not appropriate for this task.


texts: PostGIS: how to speed up casting to geography?

texts:
I have a very slow query and isolated the problem to casting my geometry type to geography. Here is my query, which basically calculates a bounding box of 15 meters (hence the need for casting to geography) around a LineString:

select st_buffer(cast(st_setsrid(lines, 4326) as geography), 15) as geometry from linesTable where id = 1

The table linesTable contains the geometries in the column lines. The cast is fine with relatively small LineStrings, but if they contain something like 500 points or more, then it takes half a second for the casting alone. The rest of the query is quick.

I need to speed this up by a factor of five or so. My concrete and to the point question is:

How can I speed up casting a complex linestring from geometry to geography? In the above query this would be the innermost part only:

cast(st_setsrid(lines, 4326) as geography

Note that st_setsrid is added for clarification only. If I drop it, then the query does not become any faster.


texts:
I've been reading a lot today about android and iphone devices being used for spatial data collection in the field and it sure seems to be that single purpose devices such as the extremely popular Garmin eTrex are a dead end. Why buy something with a tiny screen, no ability to add custom data collection interfaces, cost-plus base maps, and no picture taking when you can get all that and more for free or very cheap with a tablet/phone?

As near as I can tell so far the only real advantages a real GPS has over these other devices are ruggedness and extended battery life. Unless maybe accuracy? How does the accuracy of a GPS-enabled smart phone or tablet compare with a consumer grade "real" GPS? What about when off-grid? What about cold-start time-to-lock?


texts: How to prevent labels from one layer from overlapping features in another layer

texts: I am working ArcMap10.1 software. When I was labeling on point feature automatically label fall on line. Label must not fall on line.  I'm not aware how to solve this issues. I have attached an example image to show the problem.


texts: Not able to import LAS file into PostGIS from PDAL with filter

texts:
I'm not able to import LAS file into PostGIS from PDAL with filters.
There is no error message from PDAL, and couldn't find any relevant errors in
PostgreSQL log files, except these errors:

2019-09-19 00:00:28.804 +08 [5400] CONTEXT:  SQL function "pc_envelopegeometry" during inlining
    automatic analyze of table "postgis_25_sample.sch.table1"
2019-09-19 00:00:28.805 +08 [5400] ERROR:  function pc_envelopeasbinary(public.pcpatch) does not exist at character 27
2019-09-19 00:00:28.805 +08 [5400] HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
2019-09-19 00:00:28.805 +08 [5400] QUERY:
            SELECT ST_GeomFromEWKB(PC_EnvelopeAsBinary($1))

2019-09-19 00:00:28.805 +08 [5400] CONTEXT:  SQL function "pc_envelopegeometry" during inlining
    automatic analyze of table "postgis_25_sample.sch.table2"
2019-09-19 00:00:28.818 +08 [5400] ERROR:  function pc_envelopeasbinary(public.pcpatch) does not exist at character 27
2019-09-19 00:00:28.818 +08 [5400] HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
2019-09-19 00:00:28.818 +08 [5400] QUERY:
            SELECT ST_GeomFromEWKB(PC_EnvelopeAsBinary($1))
...

I'm using PostgreSQL version 10.10, PostGIS version 2.5.1, pointcloud_postgis ver.1.2.0, pdal 1.7.2 (git-version: Release), Windows 10 version 1903.

The LAS file is about 50GB.  I've 64GB Ram and set 300GB virtual memory.

I'm using these PDAL filters:

    {
    "type":"filters.pmf",
      "initial_distance":"0.8",
        "max_window_size":"20"
  },
 {
      "type":"filters.range",
      "limits":"Classification[2:2]"
    },
    {
      "type":"filters.divider",
      "capacity":100000000
    },

    {
      "type":"filters.chipper",
      "capacity":600
    },

What could cause the errors, how to solve the PostgreSQL errors, and where do I get the PDAL error logs?


texts: Drawing/updating feature "during" drag on OpenLayers 3?

texts:
I want to draw a line with the vertices:

initial center;
current center (on each pointer drag).

I have a blue point at the initial center and a red cross at the new center of the map on each drag.

Here is my code:

var getCenter = function() {
    return map.getView().getCenter();
};

var center = [-70000, 5500000];

var vector = new ol.layer.Vector({
    source: new ol.source.Vector({}),
    style: new ol.style.Style({
        image: new ol.style.Circle({
            fill: new ol.style.Fill({
                color: [0, 0, 255, 1]
            }),
            radius: 4
        }),
        stroke: new ol.style.Stroke({
            color: [0, 0, 0, 1],
            width: 2
        })
    })
});

var p = new ol.Feature(new ol.geom.Point(center));
vector.getSource().addFeature(p);

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM()
        }),
        vector
    ],
    view: new ol.View({
        center: center,
        zoom: 5
    })
});

var origin = getCenter();
var target = getCenter();
var coord = [origin, target];
var line = new ol.geom.LineString(coord);
var feature = new ol.Feature(line);
vector.getSource().addFeature(feature);

map.on('pointerdrag', function () {
    target = getCenter();
    coord = [origin, target];
    line.setCoordinates(coord);
});

I have created an example CodePen to illustrate what I want to achieve.

As you can see, when performing a drag, the red cross goes away from the blue point. The black line between those two only gets drawn after finishing the drag.

I want the line to be drawn during the drag.

How could I achieve this?


texts: Apply a customized function to multiple regions in Google Earth Engine Collection

texts:
I'm trying to extract the average last day of snow across a number of watersheds. The majority of the code has been provided here (https://developers.google.com/earth-engine/tutorials/community/identifying-first-day-no-snow) but I just can't figure out how to apply the customized function inside ee.image.reduceRegions for a FeatureCollection (let's call it Watersheds).
Provided code
// Define an AOI.
var aoi = ee.Geometry.Point(-94.242, 65.79).buffer(1e4);
Map.addLayer(aoi, null, 'Area of interest');

// Calculate annual mean DOY of AOI.
var annualAoiMean = annualCol.select('calDoy').map(function(img) {
  var summary = img.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: aoi,
    scale: 1e3,
    bestEffort: true,
    maxPixels: 1e14,
    tileScale: 4,
  });
  return ee.Feature(null, summary).set('year', img.get('year'));
});

With the following modifications, I would get the average value for ALL the watersheds but I want the average value for each individual watershed.
// Calculate annual mean DOY of AOI.
var annualAoiMean = annualCol.select('calDoy').map(function(img) {
  var summary = img.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: Watersheds,
    scale: 1e3,
    bestEffort: true,
    maxPixels: 1e14,
    tileScale: 4,
  });
  return ee.Feature(null, summary).set('year', img.get('year'));
});


texts: Iteration for "Selection by Location" in ArcGIS modelBuilder

texts:
I have a featuredataset of 50 Features which are supposed to have intersection with some lines and their Buffers which are in different feature classes. I have created a separate GDB for the Buffers and would like to define an Iteration in the modelbuilder which can read each feature and intersect it with all Buffers and create a new feature out of it.

I used to work with Iteration for feature classes but it is not applicable for "selection by Location".

Could you assist me how can I define such iteration for my model.


texts: How to get a collection of georeferencing control points in an ArcMap document?

texts:
I have created control points (i.e. links) using ArcMap's Georeferencing toolbar.  How can I programmatically retrieve a list of these points? I'm looking for something like an enum of the links you would see in the Link Table (which is visible through the Georeferencing toolbar).

Preferably, I would then be able to step through these points and delete any one.

I am working with the ArcObjects SDK for C# .NET.


texts: Raster calculations with values from CSV file. Iteration through attribute table

texts:
I am using GIS 3.14.16. I am working with a raster layer and a table with n multipliers in CSV format.
I would like to create n new rasters:
raster1 = raster * CSV$row1
raster2 = raster * CSV$row2
etc.

The goal is to automate this process using iteration, with n=len(csv.rows).
Is it possible to build it in Graphical Modeler?
I am attaching screenschot with the scheme of the model I would like to build.


texts:
I have four different polygons of soil data that encompass overlapping geographical areas. There is also some overlap of soil series, but each has a few unique ones:

Polygon A has Series A, B, C
Polygon B has Series A, C, D
Polygon C has Series A, B, D, E
Polygon D has Series A, E, F, G

Right now, each of them is symbolized based only on that polygon, meaning I've got different colors for the same soils in different maps.

My goal is consistent symbology for all the layers so each series always has the same color. The end product is about eight maps, and I need the same representation for all.

Since no layer encompasses the entire soil range, however, I can't copy symbology from any layer to the others without some soils being absent. To make it even more annoying, the different polygons have different naming structures for the soil types (some are abbreviated, some are broken down into additional detail [e.g. sandy clay loam vs. sandy loam], etc.).

Is there any way in ArcMap (10.1) to create a symbology from multiple layers/features at once, to cover all possible categories A though G?


texts: Disperse Markers in ArcPy or ModelBuilder?

texts:
So I have a set of around 100 clusters that are all contained in a different shapefile. I want to use the disperse markers tool on all of them, but I don't want to go in an individually use the tool on each cluster (Probably going to have to do the same process on larger sets later).

I tried using python and ModelBuilder, but the problem I keep running into is setting the representation of the layer. This is easily done through the UI by right clicking on the feature and then clicking on Convert Symbology to Representation, but this process seems much more difficult in python and ModelBuilder. I tried using variations of the representation tool (add representation, calculate representation, set layer representation).

I'm not too sure if there is a tool that can do this or some python code that may be able to do this better?


texts: Losing commmon data when joining excel and shapefile

texts:
When using QGIS 2.14.4 (Essen). I am trying to join data about North Sea Platforms with their corresponding licence block in the North Sea. My data is in excel, and contains the Block ID (see image 1), which matches the Block ID of a shapefile of the blocks themselves.

I am able to convert the excel into a csv, and then add it to QGIS Project as "No geometry (attribute only)" and then join it to the corresponding Block shapefile, allowing me to display my data by block, as desired. However, I have noticed that when two entries share the same block eg. "Alwyn A" and "Alwyn B" the resultant joined shapefile only contains information about one of them (see image 2)

Is there any way around this problem? I had thought of splitting the excel to make two shapefiles where repeats occurred, and then merging, but this seems like alot of effort, where a simpler solution may be out there.


texts:
How can I count the colorwise cell amount of a layer file? The attribute table is hidden in the menu.
Is there any tool to do that?


texts: How to change font color in one column of attribute table in Composer reports?

texts:
How do you change the font color in just one column of the attribute table included in QGIS's Composer? I would like to have the numbers listed in color red just for one column out of four columns in the attribute table.

I don't see any options to do this under Composer > Items> Item Properties> Attributes menu. (using QGIS version 2.8.2)


texts:
How to manage the symbols of the one way street a representation?
I have a field where each link is classified by one way o bidirectional.
I would like that each link read this value to insert an arrow indicating the direction of travel
My field type has this value: bidirectional - north - south - west - east (depending on the prevailing direction)
Update: I can also set a boleaan(real 0/1) field to manage the two directions, but I need to show a symbol in the middle of the link that read the orientation (0°/180° or +/-90°) from this field
In the image below a Googlemaps image with the arrows that show the direction of travel:


texts:
I'm using OpenLayers to show my map with some layers. I've a field where I can add some filters with CQL code. I can successfully add a new filter to my map, and only the right points appear in my map but when I try to remove the filters it's deleting my Layer instead of my previous CQL_Filter result, so is deleting ALL the points from the map.

//My var map is previously declared
var layerFromMap = map.getLayersByName(input.name)[0];
//If my input is empty then I should delete the filters, like resetting the search
if (input.value.length > 0) {
    layerFromMap.mergeNewParams({ 'CQL_FILTER': input.value });
    layerFromMap.redraw();
} else {
    delete layerFromMap.params.CQL_FILTER;
    layerFromMap.redraw();
}

When input is empty layerFromMapparams.CQL_FILTER is deleting my layer from the map, so all the points are deleted. Why is this happening? How can I only remove the filter so all the points appear again?


texts:
In R, I need to calculate the distance between a coordinate and all the other coordinates.
I have a dataset similar to this:

ID Morph  Sex    E   N
a    o     m    34  34
b    w     m    56  34
c    y     f    44  44

In which each "ID" represents a different animal, and E/N points represent the coordinates for the center of their home range.

I'd need to calculate the distance between each center and the centers of all the other animals. Distance with himself should not appear.

Ideally, I need to obtain two different data frames: an exhaustive one with repeated distances (a to b & b to a):

ID.1 Morph.1  Sex.1  ID.2  Sex.2  Morph.2  Dist
a       o        m    b     m      w       12
a       o        m    c     m      y       22
b       w        m    a     m      o       12
b       w        m    c     f      y       13
c       y        f    a     m      o       22
c       y        f    b     m      w       13

And another data frame with just the non repeated distances:

ID.1 Morph.1  Sex.1  ID.2  Sex.2  Morph.2  Dist
a       o        m    b     m      w       12
a       o        m    c     m      y       22
b       w        m    c     f      y       13

It is very important to me to keep the information of Sex and Morph in the new data frames! I am interested in comparing the distances among morphs (o-o, o-w, o-y, etc).

I've been doing it in Excel using Pithagoras but I am sure there must be a relatively simple code which could help me to ignite the process in R.


texts: Updating Polygon Geometries in master shapefile with detailed boundaries using ArcGIS Desktop?

texts:
I have a shapefile which contains the approximate boundaries for sites that are being worked on (thousands) in polygons, each with a reference number and additional detail in the attribution.

I have a collection of additional shapefiles with the detailed boundaries, the reference number but not accompanying information. These each have hundreds of boundaries.

What would be the best / quickest way to update the mastershapefile with the detailed boundaries, still maintaining the additional detail in the attribution contained in the master shapefile?

I could work through each shape manually, but can't get my head around how to run this as a batch.


texts:
I am using PyQGis API to create a custom application.

How do I allow users to select features with the mouse? Is there a mapTool for this? Surely with a product as powerful as QGIS there must be.


texts: QGIS 3.4.4 snapping does not lead to expected topology and st_relate result

texts:
I try to snap a line to a polygon boundary. This is working fine, due using the snapping option in QGIS (snap mode vertex and segments & tolerance 0 pixels)

Line 1 is snapping to the polygon segment.

If I export these geometries to PostGIS and do a ST_relate test with line1 and polygon, I don't get the expected result:

I would aspect that touching is true

   SELECT  st_touches (line.geom, poly.geom) FROM snapping_line line,snapping_poly poly;
   SELECT  st_relate (line.geom, poly.geom) FROM snapping_line line,snapping_poly poly;

But touches is false and st_relate is completely not what I'm expecting with snapping: 10F0FF212

I get the correct result if I snap the line to polygon vertexes (Line 2) with touching true and correct st_relate FF1F0F212

So now, where is the error? Is snapping not working correctly? Why is the result different between vertex and segment snapping. I am using QGIS 3.4.4-Madeira on Windows and Ubuntu


texts: Styling overlapping polygons of different layers

texts: I have 3 layers that sometimes have overlapping features. I need to use a thin outline style to interpret the area within the polygon and want to be able to see on the map when the features overlap. Think of it as 3 identical polygons. How can I style this?


texts:
For my class project, I need to analyze solar potential within a neighborhood at Wilhelmsburg. We need the Photovoltaic PV area, the amount of radiation taking into account the shading and the projected electricity yield of the PV system. I am using the website https://geoportal-hamburg.de/geo-online/#.  The only way was to check one by one to get the PV data. The problem is that there are 649 buildings. And doing it one by one is taking a lot of time. We can take data by using the WMS address.   But We got it as a raster file. It's good for visualization.  But not good for calculation.
It would be better if we can get this data as a vector file. Without vector file we cannot calculate. We need those data as attributes within a vector layer. Thus, you can do all sorts of calculation. Otherwise, we have to continue counting them which is too much time-consuming, because the data is very important.


texts: How to create a TIFF file using GDAL from a numpy array and specifying NoData value

texts: I have some tiff files I obtained from splitting a bigger tiff file using gdal translate. The border tiff files have some nodata values, which I can see in numpy as number 15. I want to manipulate these files as numpy arrays, and then write them back to tiff file. How can I write a numpy array to tiff file, while mantaining the number 15 as the nodata value, so they can be read propertly is GIS programs?


texts: Why are my points being flipped or mirrored when I convert from DMS to DD?

texts:
I have a collection of crash data that provides latitude and longitude in degrees, minutes, and seconds that I'm converting to decimal degrees using:

(degrees + (minutes/60.0) + (seconds/3600.0))

This appears to do a proper conversion, however the points (wgs84?) appear to be mirrored.  I feel like there's probably a simple bit of additional Math required to "fix" this, but I haven't been able to accomplish it. Any idea how I might do the conversion?

Here's a crude image illustrating what's happening:

Here's another illustration showing the points in context.  I'm assuming the shift in alignment is likely a projection issue, but the "flipping" effect is something different.

The data is provided as a CSV file, so I'm doing the conversions in Ruby.  I've used QGIS 2.2.0-Valmiera to render the images.

Here is the raw input of degrees, minutes, seconds and the output after the conversion.

45.0, 31.0, 22.9483924 = 45.52304122011111
-122.0, 40.0, 4.0819802 = -121.33219944994444

Additionally, here's the entirety of the conversion process in Ruby:

def to_decimal(row, type)
  degrees = (type == :lat ? row[:latitude_degrees] : row[:longitude_degrees]).to_f
  minutes = (type == :lat ? row[:latitude_minutes] : row[:longitude_minutes]).to_f
  seconds = (type == :lat ? row[:latitude_seconds] : row[:longitude_seconds]).to_f

  flip = type == :lon ? -1 : 1
  flip * (degrees + (minutes/60.0) + (seconds/3600.0))
end

Here's how I'm importing the CSV to QGIS (you can see some sample data in the screenshot):


texts: Removing polygons from sf simple feature collection

texts:
I'm working with a shapefile of national parks in Costa Rica, and I'm trying to get rid of the marine reserves and national parks consisting of water.

To download the dataset:

#create data driectory
    if (!dir.exists('data')) {
  dir.create('data')
}

#download dataset from ProtectedPlanet
nps_URL <- 'https://www.protectedplanet.net/downloads/WDPA_Feb2020_search_costa_rica_Costa_Rica,National_P_bf6ceea4036fb52f7e08497edb2f615999e7fffcbc1e2335dd6dbb78f3fd06a5?type=shapefile'
nps_filename <- './data/WDPA_Feb2020_search_costa_rica_Costa_Rica,National_P_bf6ceea4036fb52f7e08497edb2f615999e7fffcbc1e2335dd6dbb78f3fd06a5-shapefile-polygons.shp'
if (!file.exists(nps_filename)) {
  download.file(url = nps_URL, destfile = './data/nps.zip', method = 'auto')
  unzip('./data/nps.zip', exdir = './data')
}

#load the national parks as 'natparks'
natparks <- st_read(dsn = nps_filename)

Next, I have a vector containing the name strings of the national parks to be removed:

marinereserves <- c("Marino Ballena","Cahuita","Manuel Antonio","Tortuguero","Santa Rosa","Isla del Coco","Marino las Baulas de Guanacaste")

Then when trying to remove them, I thought of simply subsetting them by:

natparks_sub <- natparks[natparks$NAME == marinereserves,,drop=TRUE]

However, this removes nearly all national parks. Strangely enough it only keeps 2 parks, both of them are also in the marinereserves vector.

How can I properly remove the mentioned national parks from the simple feature collection?


texts:
I am using Ubuntu 16.04 and installed QGIS 2.18 following the official install details
QGIS Installers using link Release with ubuntugis dependencies
http://qgis.org/ubuntugis.
I have python 2.7.12 on my system. Every time I start QGIS, the error message is :

Couldn't load SIP module.
Python support will be disabled.
Traceback (most recent call last):
  File "", line 1, in
ImportError: No module named sip
     Python version:
        2.7.12 (default, Jul  1 2016, 15:12:24)
        [GCC 5.4.0 20160609]

        QGIS version:
        2.18.0 'Las Palmas', b7fe4bb2

        Python path:
        ['/usr/share/qgis/python', u'/home/user/.qgis2/python', u'/home/user/.qgis2/python/plugins', '/usr/share/qgis/python/plugins', '/usr/local/lib/python2.7', '/usr/local/lib/python2.7/plat-x86_64-linux-gnu', '/usr/local/lib/python2.7/lib-tk', '/usr/local/lib/python2.7/lib-old', '/usr/local/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/site-packages']

on the python IDLE i can to import SIP module.
any idea how to fix that ?

UPDATE

if i run cat /etc/apt/sources.list

show me (only related to QGIS):

.
.
.
deb http://qgis.org/ubuntugis xenial main
deb-src http://qgis.org/ubuntugis xenial main

update i have python 2.7.12 and sip version 4.17,qt version 4.8.7,pyqt version 4.11.4
after reinstall qgis now i have that error:

Couldn't load SIP module.
Python support will be disabled.
Traceback (most recent call last):
  File "", line 1, in
ImportError: /usr/local/lib/python2.7/site-packages/sip.so: undefined symbol: PyUnicodeUCS2_AsUTF8String
Python version:
2.7.12 (default, Jul  1 2016, 15:12:24)
[GCC 5.4.0 20160609]
QGIS version:
2.18.0 'Las Palmas', b7fe4bb2
Python path:
['/usr/share/qgis/python', u'/home/user/.qgis2/python', u'/home/user/.qgis2/python/plugins', '/usr/share/qgis/python/plugins', '/usr/local/lib/python2.7', '/usr/local/lib/python2.7/plat-x86_64-linux-gnu', '/usr/local/lib/python2.7/lib-tk', '/usr/local/lib/python2.7/lib-old', '/usr/local/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/site-packages']

('Qt version:', '4.8.7')
('SIP version:', '4.17')
('PyQt version:', '4.11.4')


texts: I know this is a very simple question, but how do I make a virtual raster in the latest version of QGIS (3.0)? The help files are for an older version of the program and they seemed to have some things around in the menus since they wrote them


texts: Converting DXF file to shapefile by exporting only certain fields using ArcPy?

texts:
I have a problem converting multiple DXF files into shapefiles. In particular I use the script that is shown below for conveying multiple files. But the output shapefile is very large due to the number of attributes. I would like to be able to export only the Elevation field (I am interested only to points entities) in the shp file in output. But I cannot use the FieldMapping class effectively. I have an elementary knowledge of ArcPy.

What should I do?

import arcpy
import os
from arcpy import env
env.overwriteOutput = True
in_workspace    = arcpy.GetParameterAsText(0)
outFolder       = arcpy.GetParameterAsText(1)
env.workspace   = in_workspace
datasets = arcpy.ListDatasets()
for dataset in datasets:
    filesets = arcpy.ListFiles(dataset)
    for f in filesets:
        print f
        InPoint = os.path.join(env.workspace,f,"Point")    # the CAD points
        SaniF   = f.replace(" ","_").replace("-","_").replace(".","_")  # remove invalid chars..
        if int(arcpy.GetCount_management(InPoint).getOutput(0)) == 0:
            arcpy.AddWarning("- No points")
        else:
            outPointShape = '{}_Point.shp'.format(SaniF)
            arcpy.AddMessage('- Exporting Points')
            myfields = ['Elevation']
            mapS = arcpy.FieldMappings()
            for field in myfields:
                map = arcpy.FieldMap()
                map.addInputField(inPoint, field)
                mapS.addFieldMap(map)
            arcpy.FeatureClassToFeatureClass_conversion(InPoint,outFolder,outPointShape,"#", mapS)


texts: conditional statements for raster in R - land surface temperature

texts:
Im trying to calculate Land surface temperature estimating with Landsat 8 images, but i got stuck in one place and i cant move on. I m following split window algorithm. I calculated NDVI  and I have to apply criteria to calculate land surface emmisivity. here my script in R studio :

#Ndvi
NDVI=(Palfa.B5-Palfa.B4)/(Palfa.B4+Palfa.B5)
NDVI <- NDVI@data@values
#vegetation fraction
NDVImin=0.2 #pre set value
NDVImax=0.5 #pre set value
Pv <- ((NDVI-NDVImin)/(NDVImax-NDVImin))^2
#roughness of the surface
F=0.55 #pre set value of geometric factor
C10 <- (1-es10) * ev10 * (1-Pv) * F
C11 <- (1-es11) * ev11 * (1-Pv) * F
#LSE for band 10
e10 <- if (NDVI<=0.2) {-0.047 * Palfa.B4+0.973}
       else if (NDVI>0.2 && NDVI<=0.5) {ev10 * Pv+es10 * (1-Pv)+C10}
       else (NDVI>0.5) ev10+C10

#LSE for band 11
e11 <- if (NDVI<=0.2) { -0.026 * Palfa.B5+0.984 }
       else if (NDVI>0.2 && NDVI<=0.5) { ev11 * Pv+es11 * (1-Pv)+C11 }
       else (NDVI>0.5) ev11+C11

something is wrong bcs i got this :

  Warning message: In if (NDVI <= 0.2) { :   the condition has length >
  1 and only the first element will be used ev11+C11

I also tried with ifelse but then my result is vector but i need to export raster
have u got any ideas how can i solve it?
Attached  you will find criteria :


texts: Finding centrelines from polygons in QGIS?

texts: I was wondering if there is a way to create centre lines from polygons in QGIS? I have a polygon layer, but ideally need the centrelines as a separate file.


texts: Polygon Centroid Options and Appropriateness (FeatureToPoint, Labels, etc)?

texts:
I've found a couple of other Stack Exchange questions related to this, but they were not satisfying.

For example:

How does ArcGIS Desktop calculate INSIDE point for Feature to Point?

does not tell me how the "Inside" point is determined. While the ArcObjects may refer to centroids and label points, I have at least 4 outcomes in my situation that vary differently, wildly in some cases.

I have unidentifiable points inside zip codes, so I'm assigning each one to the centroid of that zip code. Then I will calculate straight line distances from these centroids to various locations of interest. Therefore, my choice of centroid is important. Some zip codes are also of odd shapes and may possibly be discontinuous. At least that is what I assume is generating some of the odd behaviors.

Now, I have found four different ways to identify possible centrality of these polygons.

The polygon geometry contains a centroid information that can be directly extracted (and supposedly is used for labeling?).
The labeling location of the polygon (which I found does not always match where the above point is located--hence, the confusion).
Feature to point method with 'centroid' option.
Feature to point method with 'inside' option.

Some of these or most of these may overlap. In some zip codes, they are all completely different. I would have expected 1. and 3. to be similar and maybe even 1. and 2. to be identical, depending on how labeling is set up (I'm using default).

My question is not to ask "so how are these calculated?" That's neither here nor there since their algorithms are proprietary, and there's no point in trying to figure it out.

The real question is practical i.e. how should one go about determining the appropriate centroid to use?

My conclusion is the feature to point method with the 'inside' option will be the most appropriate because it is guaranteed to reside within the polygon boundary. I also assume it will be relatively central. But the descriptions about this are nebulous e.g. consider a crescent shaped polygon. The centroid would be outside the crescent but in the middle of the encapsulating envelope. Would an "inside" option literally just "push" it over into the crescent? We may not know the algorithm, but does the "inside" option actually identify centrality? Is it better than 1.? Maybe I'm missing something in my approach to this problem? Can the projection impact these comparisons? Is there a post-processing review of the polygons I should consider?


texts: No geometry when importing points from CSV

texts:
I'm having problem importing a point layer from Excel. I've saved the Excel file as a CSV and added it to QGIS (2.18) through 'Add delimited text layer' I've checked in semicolon and choosen the right coordinates.
That render a layer with correct attribute data - but no points. What am I doing wrong?


texts:
Im trying to run Focal statistics on a list of parameters, where the radius changes for each input parameter and then save the resultant DTM with the radius attached in the name. I can do this for one file, however trying to implement this within a for loop is I get the error NbrVircle object does not support item assignment. My code is below.

import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "D:/Kieran"

Parameters=[50,100,150,200,250,300,350,400]
# Set local variables
inRaster = "D:/temp/DTM.asc"
    # Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

for i in Parameters:
    neighborhood[i] = NbrCircle(i, "CELL")
    # Execute FocalStatistics
    outFocalStatistics[i] = FocalStatistics(inRaster, neighborhood, "MEAN","")
    # Save the output
    outFocalStatistics.save("D:/temp/DTM"+[i]+".asc")


texts: Full utilization of CPU-cores with gdalwarp

texts:
I need to process many large GeoTiff files with gdalwarp.
For this I have used the following command:

gdalwarp -r cubic -s_srs EPSG:25832 -t_srs EPSG:3857 -cutline shapefile.shp -crop_to_cutline -dstalpha file1.tif file2.tif /path/to/output/out.tif

The problem here is that only one CPU core is used and the calculation takes a very long time. Now I've seen that gdalwarp has the option -multi.

So I tried the following command:

gdalwarp -multi -r cubic -s_srs EPSG:25832 -t_srs EPSG:3857 -cutline shapefile.shp -crop_to_cutline -dstalpha file1.tif file2.tif /path/to/output/out.tif -wo NUM_THREADS=4

Now, gdalwarp uses 4 cores, but the utilization per core is only at 25%.

Is there a way to speed up the process or to fully utilize the cores used?


texts: Leaflet coordinate mapping wrong with OpenStreetMap when switching to EPSG 4326

texts:
I'm trying to use leaflet to add an equi-rectanglular map layer on a scatterplot of latitude and longitude points.  When I create a tile layer with OpenStreet map and option 'crs: L.CRS.EPSG4326', I get two errors: first of all the mapping transformation is wrong (for the tiles I get), and I get tile loading errors like ("GET http://b.tile.osm.org/2/7/0.png 404 (Not Found)".

Removing the crs option to go back to the default mapping fixes both of these.  Is there some additional setting or transformation that I am missing?

Code snippet:

this.mapObj = new L.Map( plotElementIDinDOM, {
            zoomControl: false,
            dragging: false,
            touchZoom: false,
            scrollWheelZoom: false,
            doubleClickZoom: false,
            tap: false,
            crs: L.CRS.EPSG4326 // L.CRS.EPSG3857 (default)
        } );
var OSM2 = L.tileLayer( 'http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
            attribution: '&copy; OpenStreetMap contributors'
        } );
OSM2.addTo( this.mapObj );
this.mapObj.fitWorld();

I am trying to map points with:

this.mapObj.latLngToContainerPoint( latLngCoords );


texts: Deleting outlines when merging two maps using ArcMap?

texts:
I downloaded two digital maps / datasets. When I downloaded the two maps, the data have  outlines. After merging two datasets, it has still outline as per the attached image.

I heard If I want to delete the outlines, I have to enter the (Open the Attribute Table -> Select by Attribute) section(?) and the delete geographic feature codes about outlines. I don't know concrete processes.
I have geographic feature codes for topographic map in Korea if that helps.

I am using ArcGIS 10.3


texts:
I have a point layer, and have got the x,y and z coordinates from 'extract values to points'. The x and y values appear fine in the attribute table but the z values are all 0. I have about 150 points, so I do not particularly want to go through each by way of info point!

Can anyone provide some advice on how to get the true z values?

I have posted a question previous to this one which is similar Extracting Z data from point layer using ArcGIS Desktop?

Slightly different as I can now retrieve and  X and Y, but Z is still a mystery.


texts: DrawPolyLine not drawing completely on the right layer

texts:
I zoom into a map, then do a screenDisp.DrawPolyline(polyline);  Now the line is in the correct place, but the map that is the background of the polyline object is actually the initial map and not the zoomed in map.

Here's a screenshot: http://i11.photobucket.com/albums/a173/iteration_x/LineBehavior.jpg

What you are seeing is a really zoomed in map of minnesota with most of the layers turned off, and that interior picture is part of the West Coast from the initial map load.


texts:
I have a French Postcode map which looks fine on screen. It's CRS is EPSG 2154 - RFG93/Lambert 93. Pretty much any French shapefile when added looks like it's leaning back suggesting to me, I'm no expert in this area, that the projection is wrong. I can see both views for each layer when I use zoom to layer on the layer in question but they are never together.

I have tried set layer and specifying the Lambert CRS and save it and open it. Still the same!

I also have a French Gazetteer text file which has lat/long. When I import that - same problem. Seems to revert to a CRS/projection/thingy that means it's leant back in the picture, and I can't see it on my main map.

Have had same issue with UK data so it's something I am doing wrong rather than the data.

Where am I going wrong?

Using QGIS Las Palmas


texts: QGIS - adding timestamps to spatialite data during data collection

texts:
I am trying to do field data collection using QGIS.  One of the features that I regard as necessary for data management is timestamps.  For example, a timestamp can indicate that an inspection record has been updated.

I have figured out how to add triggers using Spatialite to add a timestamp on update, however the timestamp is added when the edited records are saved and not when the identify dialog box is closed.  This means that timestamps can be out by many minutes.  One solution would be to force the updated record to be saved by using a custom form with an OK button (Can do the form with Qt-Designer but don't know how to add the save command to the button yet.)  However my preferred solution would be to take a timestamp for new records, when the point or line is created.  The reason for this is if I am surveying a road, I can stand on the road, collect the point and them move to the side to complete data collection.

I have been unable to work out how to get Spatialite to run a trigger that puts in the start of feature editing.  Perhaps this is also something that needs to be done by code in a custom form.  Is there an On-dirty event or signal in QGIS that indicates when someone starts editing a record and that can be used to capture the time?  Could I also ask what the Python code (and what to hook it to) that would put such a timestamp in the timestamp field.

AndrewM


texts:
I'm trying to calculate subglacial sink volume by using the rgrass7 package to utilize GRASS functions in R. However, the sinks found by the function seem to be very small:

When I run the similar functions in ArcGIS, the result makes much more sense as there are all kinds of sink sizes:

(These two images are only examples, I have verified the same results with multiple glacier rasters.)

Are there noteworthy differences between the two routines? Does GRASS not fill sinks to the pour point? Can I somehow give more parameters to r.fill.dir like the fill height etc.?

My code is:

initGRASS(gisBase = "/usr/lib64/grass76/",
          home = tempdir(),
          mapset = "PERMANENT",
          override = TRUE)

# modify current mapset with projection
execGRASS("g.proj", flags = "c", proj4 = "+proj=aea +lat_1=25 +lat_2=50 +lat_0=37 +lon_0=87 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
execGRASS("g.mapset", flags = "c", mapset = "new_mapset")

### run GRASS commands
execGRASS("r.in.gdal", input = myraster, output = paste0("GRASS_raster"),
          flags = c("overwrite"))
execGRASS("g.region", raster = "GRASS_raster")
execGRASS("r.fill.dir", input = "GRASS_raster", output = "GRASS_filled_DEM", direction = "flow_test", flags = c("overwrite")) # fill sinks
execGRASS("r.out.gdal", input = "GRASS_filled_DEM", output = paste0(output, "filled_DEM.tif"), flags = c("overwrite"))


texts: Pressing for example on listWidget that it appeared in TextEdit?

texts:
I'm doing my raster calculator in Qgis, I'm stuck in the moment how to link my listWidget with textEdit

An example like this is done in a raster calculator qgis

I want also listWidget Raster Chanel associate with expression textedit

A question as at pressing for example on 2012 that it appeared in Text edit?


texts: Downloading image/exporting image to google drive from google earth engine with color palette

texts:
When I download Analysis results from Google Earth Engine, the image is exported in to google drive but it does not have a color palette. How do i add color palette with the exported image using earth engine.

var landcover_roi = composite.clip(ft);
var ndvi =landcover_roi.normalizedDifference(['B8', 'B4']) ;
Map.addLayer(landcover_roi,{bands: ['B6', 'B4', 'B3'], min: 0, max: 0.3}, 'From Fusion Table');

// Make a palette: a list of hex strings.
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];

// Center the map
Map.setCenter(87.32,22.40, 6);

Map.addLayer(ndvi,{min: 0, max: 1, palette: palette},'Sentinel-2 NDVI')
Export.image.toDrive({
  image:ndvi,
  description: 'ndvi',
  scale: 30,
 // region: geometry,
  maxPixels:34089663741
  });


texts:
I need to access the data collected in Survey123. I reached out to ESRI and they have no training available on this topic and the documentation was lacking at best. I saw that there was a python script mentioned in the documentation. Would prefer to bypass python and interact with the REST API directly as I am trying to automate a scheduled export of the data to a reporting tool.
If the REST API can't be used are there any other ways to extract data from Survey123 that can be automated?
Based on what I have seen thus far accessing the data from the REST API for Survey123 data is poorly documented and obfuscated. I would not recommend using Survey123 if the goal is to access the data collected.


texts: oL v. 2.13: How to stop Zoom and Pan removing selected markers from array?

texts: I'm new to oL and I have an array of selected markers that is emptied each time the map redraws itself, specifically at the end of zooming or panning. I can't find where this is happening in the code and wondered if it is a built in feature of oL 2? How can I keep my selected markers selected?


texts:
I have a problem importing OpenWeatherMap to QGIS. I am new to this type of importing. I have searched old topics, but it seems that things have changed now.

Basically, I have a map in QGIS, and I want to overlay the existing map with a real-time weather data.


texts: GNAF query for all points in Turf denoted by list of points?

texts:
I'm using Australia's  gnaf geospatial info which is loaded into postgres with the postgis extension. It is a SQL database . Gnaf has a table of address with long,lat and i have a list of points like , (-1,1) (1,1) (1,0) .... n times .
I need to query for all the points that fall inside the polygon formed by the points.

Can anyone suggest an algorithm that can be used to build the query.?

An example query to extract all points while applying restrictions

SELECT * FROM Gnaf WHERE longitude > x AND latitude < y


texts: Modifying styling on infowindows (Carto)?

texts:
I'm working on a map with lots of infowindows, many of which are too tall for the embed window, which requires some awkward scrolling in order to exit out of any given card. Is it possible to modify the styling of the infowindow on Carto to make the pop-up cards horizontal? Or to change the dimensions of the map such that it is square rather than widescreen?

https://rszypko.carto.com/builder/050f7b58-0dae-11e7-86df-0ee66e2c9693/embed


texts: Why does ArcMap connect to a Personal Geodatabase (PGDB) even when no feature classes (or tables?) are referenced in the MXD?

texts:
I have an odd question from a user in the same goverment that I work for.

When he opens a .MXD-project in ArcMAP a Personal Geodatabase also opens that is not part of the project and does not have any connection join/relates to the project (an LDB-file is created an points to his computer). He removed all his OLE DB-connections and saved the project under a different name but the problem still occurs.

Does anyone have a clue of what's going on here?

I myself think the user opened the PGDB in Access and did something the he shouldn't have.


texts: Why does intersect not calculate polygon areas?

texts:
this is strange: I intersected 2 shape-files and the resulting shape does not contain the Shape_Area and Shape_Length of the intersected polygons. Instead the field Shape_Area contains the polygon-area of one of the 2 original Layers.

Recalc of the geometries with vector->geometry->add/export column seems to add the areas. That is quite a drag and very slow for my 350.000 polygons. So I am hoping for a systematic approach that adds the areas while generating the intersect.

Is intersect done by QGIS (qgis.bin) or by a gdal-tool? Either way: is there a "calculate areas"-option hidden somewhere?

THX, stn


texts: How can I get a .asc file to load into QGIS?

texts:
I got some files from an ag plane that are in ASCII format. I cannot seem to be able to get them to load into QGIS without the

This is a link to the file I am trying to use.

https://drive.google.com/file/d/0B7IZBVWB4R68Q1dBNlJaaUFUeDQ/view?usp=sharing

I was wondering if it was a problem with the file itself or is there an easy (or hard) solution to getting QGIS to read the file correctly?


texts: Packing squares polygons in a raster grid

texts:
In the green area of the raster I have to fit as many square polygons as possible, of three different sizes:

3x3 cell size
2x2 cell size
1x1 cell size

With the fact that the approach should be greedy, so first fitting the biggest polygons, then the middle sized, and then if there is any space left, the smallest size. It does not have to be an exact approach but at least an approximation.

Considering this is just a test area for the algorithm, which should eventually perform this search on bigger real case rasters, where basically we have objects of three different sizes and are trying to place them in the allowed zones of the raster.

Having said that, I am aware it's not at all an easy task. I am using Qgis with R scripting to try to formulate this search.

EDIT: By moving a 3x3 window square and suming up the cell values, it can be determined where are the possible locations for the squares:

So, wherever there is a 9 means the square can fit, but the squares cannot overlap so now the question is how to find the best combination of 9 cells that do not overlap


texts: cache in geowebcache does not update when feature delete or added

texts:
I have Layer with geometry point when I added new point or delete feature point tile in geowebcache does not update if I want to saw updated  tile I must delete folder layer in geoweb cache?
what can I do?
I think that  geowebCahe automatic update when feature modify and tile is updated .is it true?
How can update tile manually in geowebcache?


texts: In ArcMap 10.4 I'm trying to manually join two nearby points on a one-by one basis by visually inspecting a single point in one layer, and associating it somehow with a point in another layer.  Using nearby points is troublesome, because I need to use subjective criteria to make that join.  For example, I have a geolocated electricty service account point, and I need to associate it with a gas service account.  They are not always the two closest points on a map, so I must use my judgement to join them together.  The last step would be to join those points to a building footprint which is a geolocated polygon.  I'm not savvy in python, so a way to click one point and then click a second point and join them would be something I initially think would be a user-friendly solution?


texts: How to use wildcard when appending PDFs using ArcPy?

texts:
I'm working on a script that (as written) produces a large number of individual pdf pages (>100), instead of hardcoding the names of each of the outputs, I was thinking of using a wildcard. However, neither of the following seem to work:

 pdfDoc.appendPages("C:\\temp\\*.pdf")

or

 pdfDoc.appendPages("C:\\temp\\"*".pdf")

Below is the full script:

mxd = arcpy.mapping.MapDocument("C:\\temp\\Graffiti_DM3.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
graffitilyr=arcpy.mapping.ListLayers(mxd, "*", df)[0]

newquery = '"Incident_Z" =' + "'10003'"
graffitilyr.definitionQuery=newquery
features = arcpy.SearchCursor(graffitilyr)
count = 1

for feature in features:
        df.rotation = "0"
        arcpy.SelectLayerByAttribute_management(graffitilyr, "NEW_SELECTION", "FID=" + str(count))
        #df.extent = graffitilyr.getSelectedExtent(True)
        df.zoomToSelectedFeatures()
        df.scale = df.scale * .001
        arcpy.RefreshActiveView()
        df.panToExtent(graffitilyr.getSelectedExtent())
        #df.scale = "1000"
        outpdf="C:\\temp\\" + str(count) + ".pdf"
        arcpy.mapping.ExportToPDF(mxd,outpdf)
        count = count+1

pdfpath="C:\\temp\\pdfbook.pdf"
pdfDoc = arcpy.mapping.PDFDocumentCreate(pdfpath)
for pdf in pdfs:
        pdfDoc.appendPages("C:\\temp\\*.pdf")
pdfDoc.saveAndClose()
del features
del pdfDoc

Any help would be much appreciated


texts: ArcGIS Pro field calculator if - else - string error

texts:
I try to calculate some fields with the field calculator using a simple if-else condition. I would like to write the values in the field volume, based on some calculations including values of fields wert and shape_area (no coordinates).

The code works, when I use just one integer/no math functions.
If I try it like it should be, I get syntax error .

The example:

import math
    def muh(bedingung, wert, Shape_Area):
       if (bedingung == "hallo"):
        value = (4.3 * math.log(!wert!)) * !Shape_Area!
       elif (bedingung == "echo"):
        value = (2.7 * math.log(!wert!) * !Shape_Area!

       return value

And expression volume = "muh(!bedingung!, !wert!, !Shape_Area!)

I can't seem to find the fault. The problem lies between numeric values and "string" (?), I guess, I have to escape or convert some parts of the code?


texts: Dynamic legendlabel for a rule-based style

texts:
I have a rule-based style that is linked to a print-composer.  In that composer I activated the atlas-function with a filter.
Because of that filter the outcome of the rule-based style varies.

Is there a way to make the labels in the legend also dynamic, so the legend shows the exact results without putting them in manually?
I'm using QGIS 2.18.18.


texts: Clipping TileLayer with vector polygon (clipping mask)

texts:
I would like to add a mask to clip layer (Tile) in a way that only the part bounded by a ol.Feature (polygon or multipolygon or layer) is visible and the rest is hidden (left transparent). The geometry is generated dynamically so I need a way to do that programatically.

I presented the question Clipping TileLayer with georeferenced polygon (clipping mask) (see).
Clipping map with vector layer (only layers zIndex == 0):

var style = new ol.style.Style({
    fill: new ol.style.Fill({
    color: 'black'
  })
});

var clipLayer = new ol.layer.Image({
  source: new ol.source.ImageVector({
    source: new ol.source.Vector({
      url: '//openlayers.org/en/v4.6.5/examples/data/geojson/countries.geojson',
      format: new ol.format.GeoJSON()
    }),
    style: style
  })
});
clipLayer.on('postcompose', function(e) {
  e.context.globalCompositeOperation = 'source-over';
});
clipLayer.on('precompose', function(e) {
  e.context.globalCompositeOperation = 'destination-in';
});

Example result: OSM (BaseLayer), TileLayer (Stamen) and vector layer

Two layers were clipped (!) - OSM and TileLayer.

I need clipping TileLayer with vector layer (support zIndex).
Example: OSM (BaseLayer), TileLayer and vector layer

One layer were clipped - TileLayer.

Does anybody have some clues or examples to achieve this?


texts: Loading GML data into Postgres with ogr2ogr

texts:
I'm working with a new OS dataset, Highways Network, trying to get it loaded into Postgres and set up in some kind of usable state.

I've loaded all the data already, but it seems some elements of what is contained in the GML files have not been included.

Here is a section of the GML:

<ram:Maintenance gml:id="id_1770MA01526959">
  <gml:identifier codeSpace="http://inspire.jrc.ec.europa.eu/ids">http://data.os.uk/id/1770MA01526959</gml:identifier>
  <net:networkRef>
    <net:NetworkReference>
      <net:element xlink:href="#usrn44302253"/>
    </net:NetworkReference>
  </net:networkRef>
  <net:inspireId>

The element that has not loaded is: xlink:href="#usrn44302253"

If I run ogrinfo on the GML, it does not even pick up the #usrn:

OGRFeature(Maintenance):131433
gml_id (String) = id_1050MA01364874
identifier (String) = http://data.os.uk/id/1050MA01364874
localId (String) = 1050MA01364874
namespace (String) = http://data.os.uk/
beginLifespanVersion (String) = 2016-03-15T00:00:00.000
maintenanceResponsibility (String) = Maintainable At Public Expense
maintenanceAuthority|ResponsibleAuthority|identifier (Integer) = 1050
authorityName (String) = Derbyshire
partialReference (String) = false
highwayAuthority|ResponsibleAuthority|identifier (Integer) = 1050
highwayAuthority|ResponsibleAuthority|authorityName (String) = Derbyshire
reasonForChange (String) = New
locationDescription (String) = (null)

Everything else listed above has been loaded successfully into postgres.

Can anyone tell me why this is, and how I can persuade ogr2ogr to load the #usrn as well as everything else?


texts: GeoServer vector tiles (.pbf) wont display due to CORS issue

texts:
I am using Mapbox Gl JS to get mapbox vector tiles (pbf) from GeoServer 2.15, which I am running with Tomcat 8.5 on a Windows Server. The Mapbox GL JS javascript is in a folder in the Tomcat webapps directory. Tomcat has been configured for CORS by editing the web.xml to include:

  <filter>
  <filter-name>CorsFilter</filter-name>
  <filter-class>org.apache.catalina.filters.CorsFilter</filter-class>
  <init-param>
    <param-name>cors.allowed.origins</param-name>
    <param-value>*</param-value>
  </init-param>
  <init-param>
    <param-name>cors.allowed.methods</param-name>
    <param-value>GET,POST,HEAD,OPTIONS,PUT</param-value>
  </init-param>
  <init-param>
    <param-name>cors.allowed.headers</param-name>
    <param-value>Content-Type,X-Requested-With,accept,Origin,Access-Control-Request-Method,Access-Control-Request-Headers,Authorization</param-value>
  </init-param>
  <init-param>
    <param-name>cors.exposed.headers</param-name>
    <param-value>Access-Control-Allow-Origin,Access-Control-Allow-Credentials</param-value>
  </init-param>
</filter>
<filter-mapping>
  <filter-name>CorsFilter</filter-name>
  <url-pattern>/*</url-pattern>
</filter-mapping>

I have also added URL rewrite within IIS which links port 80 to Tomcat/GeoServer on port 8080. When opening the map in browser on another machine everything but the vector tiles from GeoServer displays. This is the error message I get:

  Cross-Origin Request Blocked: The Same Origin Policy disallows reading
  the remote resource at
  http://localhost:8080/geoserver/gwc/service/tms/1.0.0/webmap:target_status@EPSG%3A900913@pbf/6/53/25.pbf.
  (Reason: CORS request did not succeed).

I have also noticed that X-Frame is set to SAME ORIGIN for the vector tiles.

Is there anything else missing or wrong with my setup? I would like to be able to open the web map through the URL and for the layers to be displaying? Is the CORS issue coming from IIS or Tomcat?

Edit: Just wanted to clarify that GeoServer itself is working as I can access it from another machine, through the URL. It is just the vector tiles which are generating the CORS error when viewing the map.


texts: Write array to GeoTIFF in polar projection - Python

texts:
I have done some calculations on an subset of an hdf5 file in python. I now want to write my results to a GeoTIFF to analayze in QGIS.
I have the 664 rows by 632 cols data array that I've calculated from the hdf5 subset. The HDF5 also has lat and lon subsets with the same dimensions as the data.
The file is projected in WGS 84/NSIDC Sea Ice Polar Stereographic South (EPSG 3976).
When I open the lat and lon in Python the lat and lon files are the same as each other and have values values in the 2000s. Top left corner is 2471, top right 2685, bottom right 2532, bottom left 2505. When the same lat lon subsets are opened in QGIS the pixel values have lat/lons that make sense for a map projected with the south pole at the centre.
The following site's bottom table lists the corner coordinates in lat/lon and a grid. https://nsidc.org/data/polar-stereo/ps_grids.html.
I am trying to write my data to a TIFF projected in EPSG 3976 using the answer by EddyTheB herebut cannot get it to work. I get a GeoTIFF with the data in a grid, but it appears in a nonsense area of the map in QGIS compared to other data I have loaded and correctly projected. I have tried various combinations of xmin, xmax etc using lats from the website's grid, the lat lon files in python etc but none work. I expect it is something to do with it being different with a polar projecition?
(th89 is my data array)
import numpy as np
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr
import matplotlib.pylab as plt

xmin,ymin,xmax,ymax = [lon2.min(),lat2.min(),lon2.max(),lat2.max()]
nrows,ncols = np.shape(th89)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)

geotransform=(xmin,xres,0,ymax,0, -yres) # (ulx, resx, 0, uly, 0, -resy)
# That's (top left x, w-e pixel resolution, rotation (0 if North is up),
#         top left y, rotation (0 if North is up), n-s pixel resolution)

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols, nrows, 1 ,gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  # Specify its coordinates
srs = osr.SpatialReference()                 # Establish its coordinate encoding
srs.ImportFromEPSG(3976)                     # NSIDC is 3976
                                             #
                                             #
output_raster.SetProjection( srs.ExportToWkt() )   # Exports the coordinate system
                                                   # to the file
output_raster.GetRasterBand(1).WriteArray(th89)   # Writes my array to the raster

output_raster = None

Geolocation info as follows (corners clockwise from top left):


texts: How to get an GeoExt LegendPanel with image/png instead of image/gif?

texts:
I am using a GeoExt.LegendPanel to show legends for my various layers in OpenLayers. The GeoExt.WMSLegend supports configuring which image format to use, but the LegendPanel does not utilize this causing the WMSLegend to fallback to image/gif.

So far I have manually edited the GeoExt.js minified file to make the WMSLegend use image/png as default instead of image/gif, but I would prefer to parameterize this.

Am I missing some way to configure my LegendPanel to use image/png as my WMSLegend image format?


texts:
I've got a shapefile that was made by the mother of all Unions. Because of that, it's got about 400 columns, and I want to sum about 100 of them (specifically, the "Score" from each). Is there an easier way to do it than to manually type a formula adding each score value?

Before this project is done, there will be about 600 scores that need to be added to one thing or another so any kind of automation would be helpful. I've thought about setting up some kind of iteration through the column names (they're sequentially numbered) but I don't know how to do that. I've also thought about turning off all the fields other than the scores (no fun) and then summing all values for each row (which I also don't know how to do).


texts:
I have an ESRI shapefile of the contours of my study area, and I want to convert it into a format that is acceptable by GMT, so I can make an elevation map of my study area. Is that even possible?
Would it be easier, if I had a dem file, to convert it into gmt format?

Edit: I made a .dem file. How can I convert it to GMT acceptable type?


texts: How to obtain a network completely connected in QGIS?

texts:
I have a vector of lines, in which some of them are isolated from the rest of the lines, but I need that every line is connected to at least another one, so as to obtain a unique connected path.

This is a part of the network, where there are isolated (black) lines:

I would like to know if there is a tool that could connect automatically those lines in order to obtain a unique path (for example, as I shown with the evident handmade red lines)


texts: Why aren't my hyperlinks saving in ArcGIS 10.2

texts:
I have a polygon layer and I have a joined table with a link to an excel document for each polygon. Under the "Display" tab in layer properties I check the "Support Hyperlinks using field:" box and select the joined field that my document link is in and choose "Document" as my link type, apply and OK. It works fine, the excel document opens when I click the polygon with the hyperlink tool as it should. Now when I save the .mxd file,close it and then open it again, my changes to the Hyperlinks are not saved. Any other changes I make are saving, my joined table is still there, but the "Support Hyperlinks using field:" is no longer checked.

Why would it not be saving this?


texts: How to stack multiple graphics/marks/backgrounds on a single TextSymbolizer in Geoserver?

texts:
I tried simply adding two marks to a graphic, but only 1 got rendered. I tried adding two graphics to a TextSymbolizer, but only 1 got rendered.

I tried adding two TextSymbolizers to the same rule, and it partially worked; they both got rendered, but they got rendered in different positions or stacked in the wrong order.

To fix the wrong order, I tried changing their priority, but this caused them to never render on the same position. That's probably because of the space-around VendorOption, so I set that to -1000.

This made it work in about 75% of the cases. I got them both to render, and in the correct order, but sometimes, I get one of the layers rendering in a different spot, or two different TextSymbolizers will overlap because of the negative space-around. The Text has varying length, so I can't use fixed sizes etc.

My desired result can be seen in the image below; A mark (of type square) with a stroke stacked ontop of a bigger mark (also square). I tried creating a new mark (by implementing MarkFactory etc), but I only want the stroke on the inner mark, and I couldn't find any way to make an edge "unstrokeable".

Any way to accomplish this in an easier way? Perhaps a way to get multiple strokes on a mark?


texts:
We have set up a little form with a few attributes and a photo.
Once the photo is saved, it is not possible to reopen it and to zoom in the picture.
Any tricks?


texts: QuickOSM error occurred in executing Python Code

texts:
So I'm trying to download a dataset using QuickOSM, and I keep getting the same error over and over again. I thought at first the selection size was too big, but making it smaller didn't fix it. I wonder if anyone else has come across this error and how to fix it.
Here is the dialogue that I get when I try to execute it:
An error has occurred while executing Python code:

QuickOSM.core.exceptions.NetWorkErrorException: ('Nominatim API : Download failed: Network access is disabled.', None)
Traceback (most recent call last):
  File "C:/Users/sbrow/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QuickOSM\core\api\nominatim.py", line 54, in error
    raise NetWorkErrorException('Nominatim API', ', '.join(messages))
QuickOSM.core.exceptions.NetWorkErrorException: ('Nominatim API : Download failed: Network access is disabled.', None)

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.14.16-Pi Pi, df27394552

Python Path:
C:/OSGEO4~1/apps/qgis/./python
C:/Users/sbrow/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/sbrow/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/OSGEO4~1/apps/qgis/./python/plugins
C:\OSGEO4~1\apps\Python37
C:\OSGEO4~1\apps\Python37\Scripts
C:\OSGeo4W64\bin\python37.zip
C:\OSGEO4~1\apps\Python37\DLLs
C:\OSGEO4~1\apps\Python37\lib
C:\OSGeo4W64\bin
C:\OSGEO4~1\apps\Python37\lib\site-packages
C:\OSGEO4~1\apps\Python37\lib\site-packages\win32
C:\OSGEO4~1\apps\Python37\lib\site-packages\win32\lib
C:\OSGEO4~1\apps\Python37\lib\site-packages\Pythonwin
C:/Users/sbrow/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/sbrow/New folder


texts: Copying feature from one feature class to another without creating new feature class using ArcGIS Desktop

texts:
I am using ArcGIS Desktop.
How do I copy a line feature from a line feature class to another existing feature class without creating a new feature class?
The Copy Features tool creates the copied feature in a new feature class but what I want is to add the feature that will be copied to an existing feature class.


texts: Downloading JAXA/GPM_L3/GSMaP/v6/reanalysis" hourly product in Google Earth Engine

texts:
I want to download "JAXA/GPM_L3/GSMaP/v6/reanalysis" hourly product for a relatively large area approximately 300000 square km and for a long period of time, say 15 years, since this product is hourly, the resulting tiff files will be too many and after 2 hours of waiting for the product to be ready for download from google drive the running process stops for an Internal Error message.
Do I need to aggregate hourly data to at least daily  summation to make it easier to download?
The link to my code for downloading data in stacked form is
https://code.earthengine.google.com/17a92bae24f7445113495668883d03d8. The code is also mentioned in the next lines.
var GSMap = ee.ImageCollection("JAXA/GPM_L3/GSMaP/v6/reanalysis")
.filterBounds(geometry)
.filterDate("2001-01-01", "2014-03-01");

// Choosing bands.
var Preci = GSMap.select("hourlyPrecipRateGC");

// Stacking bands.
var Preci_stack = Preci.toBands();

// Converting band names which is a list of dates to feature collection.
var bandNames = Preci_stack.bandNames();
var fc = ee.FeatureCollection(bandNames.map(function(point) {
  return ee.Feature(null, {'value': point});
}));

// Exporting the generated feature collection to csv.
Export.table.toDrive(fc, 'GSMap_Hourly_Precipitation_test');

// Exporting generated images to google drive.
Export.image.toDrive({
  image : Preci_stack,
  description : "GSMap_Hourly_Precipitation_test",
  scale : 10000,
  region : geometry,
  crs : "EPSG:4326",
});```


texts: Checking if an attribute matches a lookup table

texts:
I have a point file representing plants (Species 1, Species 2, etc.) and a polygon file of ecozones (Zone A, Zone B, etc.). I also have a lookup table for which zones each plant species is known to grow in, which looks like:

Species A, Zone 1
Species A, Zone 2
Species A, Zone 3
Species B, Zone 2
Species B, Zone 3
Species B, Zone 4
etc.

My goal is to determine whether or not each plant point is inside a zone it is known to grow in. By using a spatial join, I am able to apply the zone code to each point, but how can I relate this value to the lookup table to determine whether the zone assigned is in the lookup table? Is there something I should be doing before the spatial join?


texts: Tiling/Caching Image Services with ArcGIS Server

texts:
The company I work for uses online map services while working with our clients.  We publish map documents as well as image services and have had no issues until recent.

We are able to publish/create a service for map documents containing vector data only with no issues but are now having issues with our raster services.  We are having two seperate problems:

1 - We have geodatabases with raster mosaics of LiDAR data which we can successfully publish however they will not tile/cache; and

2 - We have a geodatabase with a raster mosaic of imagery which we can't publish or tile.  When publishing we get the error "Packaging Succeeded, Publishing Failed"

I've tried to do run both situations from various computers in the office as well as directly off the server and nothing seems to work.  This issue only started a few weeks ago.

Any suggestions on what might be the issue?


texts:
I want to do the NDVI with images landsat 8 in R, radiometric correction, radiance and reflectance, I make it from the models of the USGS (model here)

But taking the image (ND) for the calculation of the radiance and this for the calculation of the reflectance, it gives me errors that I then show pixel values within the image in a cloudless zone.

Evidently the values of the radiance are wrong the calculation of the NDVI gives wrong. The research has led me to know that it is something related to the 16 bits in which the image is. These are the pinxel values of the original image.

As I can solve this problem in R and that the radiance d values give me good to be able to do the NDVI and me between -1 and +1...

Something similar was treated in another post but was solved given that ENVI internally corrects these values without doing intermediate processes (Ref:Processing landsat 8 in ENVI)

values ND, Radiance, Reflectance


texts:
i need to print a detailed map, on a street level zoom of Metro Manila. i figured to do this i need to download map tiles, but that will be a breach on map usage in OSM. then i found this website, http://metro.teczno.com/#manila, they provide pre-downloaded map extracts of certain areas, with Metro Manila being one, just what i need.

Now what i have to do is get a hold of some GIS application, and that will be QGIS, I've done some research and it turns out QGIS offers a plugin that prints maps, print composer to be exact, but as i continue reading, though i haven't really tried it, print composer only prints in a certain zoom level? or paper/print size?

I need a map that is on a street level, so if im gonna use one paper, it has to be huge to accommodate the street zoom level, so im thinking if print composer provides or has a feature to cut a map on a street level zoom to multiple A4 size papers or any smaller size papers, so that i can just later on fit them all together?


texts: Does anyone know how to change the properties of a vertex to keep 3 places after the decimal place? My problem is when I go to draw a line and snap to a point, the properties of the line I am drawing are rounding to the nearest whole number instead of keeping the correct number of decimal places and the vertex positions arent staying where I need them to stay. Is there some way to fix this in some kind of settings or options menu. In the Editor Menu - Options, I have Display measurements using 3 decimal places but when I start drawing a line, it rounds to a whole number for each vertex so every vertex is just a little off of where it needs to be.


texts: How can I calculate the number of raster pixels within polygons?

texts: I have two layers, flood area (flood.img) as a raster and flow direction as polygons. Would you please help me in calculating the number of pixels in each flow direction polygon?  I want to get the percentage of flood area in each flow direction polygon. I use ArcGIS 10.


texts:
I use QGIS 3.8.0 and I'm looking for tools to assign a polygon parent depending on the value. I attach the picture, where there are several polygons with different names.

I would like to get the name of the parent (polygon with neighboring value and containing this polygon) in the attribute table.

I've tried to use aggregate function but my skills are to weak.

I've used this command:

 aggregate(layer:='layer',
 aggregate:='concatenate',
 expression:="Name",
 concatenator:=',',
 filter:= contains( $geometry, geometry(@parent) )
)

It should look like this:

Anyone have any idea how can I get this information ?


texts:
I'm in the process of adding automatic CRS conversion to a Spark application and running into some problems with the implementation.

Each row in a dataframe that is being converted will contain columns Latitude, Longitude, and EPSGCode and will result in new columns for the converted Lat/Long in WGS84.

Spark itself isn't the issue (included for context only), as wiring this up to run within a UDF is fairly simple.

Right now, I'm running into issues with a common conversion, given EPSGCode 4267 and targeting 4326, I either get an error or just get back the original coordinates.

I've tried two libraries...

OrbisGIS CTS: https://github.com/orbisgis/cts (throws error)

Proj4j: https://github.com/locationtech/proj4j (returns same points)

Here's a sample of my scala code using CTS:

class DatumConversionSpec extends FunSpec {

  val cRSFactory = new CRSFactory()
  val registryManager: RegistryManager = cRSFactory.getRegistryManager
  registryManager.addRegistry(new EPSGRegistry)

  def getCRSFromEpsgCode(code: Int): GeodeticCRS = cRSFactory.getCRS(s"EPSG:$code").asInstanceOf[GeodeticCRS]

  it("handles NAD27 to WGS84 conversion") {

    val source = getCRSFromEpsgCode(4267)
    val target = getCRSFromEpsgCode(4326)

    println(source)
    println(target)

    // throws exception here!!!
    val ops = CoordinateOperationFactory.createCoordinateOperations(source, target)

    val op = CoordinateOperationFactory.getMostPrecise(ops)

    .....

which throws the exception java.lang.IllegalArgumentException: No transformation found from EPSG:4267 to EPSG:4326.

When using proj4j I just got back the original coordinates.

Is there another library that could recommended for this or a different approach I should be taking?

note: our Spark infrastructure runs Scala 2.11.8 currently against Java 8. Upgrading to a newer Java version to use another library is not an option right now


texts: Extract proportion of raster by several rectangles

texts: I have quite big raster and I 7 restangle polygons. What I need is to extract the raster so the result is 7 small TIFF tiles ideally using QGIS (possibly also ArcGIS but it's really slow here). Clip feature doesn't work properly as the result is one big black polygon.


texts:
I have 2 sets of data:

A network data and I want to validate it with another
B network geometry which is a moving van trace data

Since the trace data is very raw without any attribution, is there a way to conclude that I have captured Anetwork wherever BNetwork is present?
I tried simple select by location, but that doesn't work. Plus the data is very heavy (1M + records).


texts: Is there a vector format for depicting continuous 2D polygons with attributes?

texts:
I'm working within the QGIS environment and looking for a 2D continuous vector attributed editable file format to work with.

My goal is to find a vector layer that has different features with no topological gaps nor overlaps among them. Is there a well-known file format that does that?

This layer needs to have its geometry easily edited in respect with its topological constraints. I prefer no post-editing processes like ArcGIS' topological fixes. What I had in mind is something like the ArcGIS' Parcel Fabric type but with a shorter and easier editing process.

----Places I've looked in----

OGR Vector Formats page.
OGC® Standards and Supporting Documents
Dr. Google. With various keywords...


texts: TableToTable tool fails in ArcObjects depending on output location?

texts:
The code I am writing for a button (ArcObjects 9.2, Visual Studio 2005) takes the attributes for selected points in a shapefile and writes them to a table:

Dim pGp As GeoProcessor = New GeoProcessor
pGp.OverwriteOutput = True
Dim pParams As IVariantArray = New VarArrayClass
...
Dim strFolderName As String = "C:\temp"
Dim strTableName As String = "TandE_data2.dbf"
With pParams
    .RemoveAll()
    .Add(pTELayer) 'A point SHP with at least one feature selected
    .Add(strFolderName) '("in_memory")
    .Add(strTableName) '("TandE_data2")
End With
pGp.Execute("TableToTable_conversion", pParams, Nothing)

As written, the process will fail at the pGp.Execute statement on the third execution, with the message

A first chance exception of type 'System.Runtime.InteropServices.COMException' occurred in SDtools.dll

appearing in Visual Studio's output window. If the output is directed to an in_memory workspace table instead of a DBF (using the values seen in the code comments), the code always runs to completion.

What changes do I need to make to the code for it to output to a DBF consistantly?