texts: QGIS3: Automating making figures for all raster layers in a group

texts:
I'm using QGIS 3.14 and I have a group ("Constant Climate" in the example") containing several raster layers, and I have created a Layout with the Print Composer.
Now I would like to make a map for each of the rasters in the group, where only the selected raster is shown, and and all the other rasters in the group are invisible. This is akin to the Atlas where the group is the equivalent to the vector layer selected in "Coverage layer", and the rasters in the group are the equivalent to the features in the selected vector layer.
Turning layers on and off by hand and then exporting the figures is cumbersome, in particular since this is a task I need to do frequently, and in many of my QGIS projects.
Is there an elegant way to do this, or does it involve writing python code (I use python a lot outside of QGIS, but I have limited experience with pyqgis)?


texts:
This arises from my own question How to handle coordinates accuracy in ArcGIS where I have tried to use the documentation entitled Using geometry objects with geoprocessing tools as a reference.

I have a table with coordinates in degrees:

I created event table and added it to the view with coordinates system 'GCS_NZGD_2000
WKID: 4167 Authority: EPSG'. I converted this single point to shapefile, defined it projection and computed coordinates of the point using 'Add Geometry Attributes' tool. This is resulting table with numbers as expected:

To replicate this in arcpy I've used this code:

corners =[[174.73,-36.76]]
p=[arcpy.PointGeometry(arcpy.Point(*coords)) for coords in corners]
arcpy.CopyFeatures_management(p, "d:/rubbish/points.shp")

I added output 'points.shp' to the view, defined projection and computed coordinates of the point using 'Add Geometry Attributes' tool. This is resulting table:

As one can see from the picture below the distance between 2 supposedly identical points is close to 10 meters:

However when I updated existing dataset with defined projection using

infc =r'd:\scratch\from_xy.shp'
outRows = arcpy.da.InsertCursor(infc,("SHAPE@","X"))
feat = (arcpy.Point(174.73,-36.76),0)
outRows.insertRow(feat)

It worked. Lessons:

Don't use examples similar to the documentation entitled Using geometry objects with geoprocessing tools
Define projection of dataset prior to any games with geometries.

Does the documentation entitled Using geometry objects with geoprocessing tools need to be revised?


texts: Wildcards in feature selection using geopandas?

texts:
I need to select the rows that are according to the following line:

Like this is for numbers:

data[data['field']<4]

I need this: where the words do are inside the name of a given row

data[data['field'] LIKE '%do%']

and this: where it starts with 01

data[data['field'] LIKE '01%']

I think it needs regex but I don't know how to do it


texts: Problems georeferencing an image in QGIS with subprocess.popen() and gdal_translate

texts:
I want to georeference a raster using GDAL in QGIS 3.10's Python console. When I use the plugin I get the following script (leave only two GCPs to make it easier to read):

gdal_translate -of GTiff -gcp 718098 -6.35036e + 06 718099 6.35036e + 06 -gcp 712827 -6.3758e + 06 712833 6.37585e + 06 "input.tif" "temporal_output.tif"

gdalwarp -r near -order 2 -co COMPRESS = NONE -t_srs EPSG: 32720 "temporal_output.tif" "final_output.tif"

These are the GPCs that I am using:

The code I'm running to try to replicate the results is as follows:

gcp = open('input.points', 'rb')
hdr = gcp.readline()

command = ['gdal_translate', '-of', 'GTiff']

for line in gcp:
    x, y, col, row, e = line.decode().split(",")
    command.append("-gcp")
    command.append(str(round(float(col))))
    command.append(str(round(float(row))))
    command.append(str(round(float(x))))
    command.append(str(round(float(y))))

command.append("input.tif")
command.append("temporal_output.tif")

proc = subprocess.Popen(command)

I do something similar with the gdalwarp command, but I think I have the problem with gdal_translate. The command creates a file "temporal_output.tif", but when I compare the geographic parameters I have:

g = gdal.Open("input.tif", 1)
gt = g.GetGeoTransform()
print(gt)
(699129.0, 2.5, 0.0, 6392937.5, 0.0, -2.5)

g = gdal.Open("temporal_output.tif", 1)
gt = g.GetGeoTransform()
print(gt)
(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)

El objecto 'command' queda de la siguiente forma (reduced the GPCs for a better reading):

print(command)
['gdal_translate', '-of', 'GTiff', '-gcp', '718098', '6350356', '718099', '6350357', '-gcp', '712827', '6375795', '712833', '6375855', 'input.tif', 'temporal_output.tif']

(I also replaced command.append(str(round(float(row)))) with command.append(str(round(-1*float(row)))), to make it 'more similar' to the plugin script, but it doesn't modify the result at all)

I'm new to python and gis, and I can't think what the code error would be.

Can anyone find the error?

Thanks in advance


texts:
I tried to add more than three rules of CSS styles, whenever I tried to add fourth rule, I can't submit the style in Geoserver 2.8.1.

Please find the CSS styles I used below.

/* Default Style less than a scale*/
*[@scale < 50000] {
   fill: #FFD9B3;
   fill-opacity: 50%;
   stroke: #D9D8D8;
   stroke-width: 0.5;
}

/* Default Style 2 */
[@scale >= 50000] [@scale < 250000]{
   fill: #FFCC99;
   fill-opacity: 50%;
   stroke: #E2E0DD;
   stroke-width: 0.5;
}

/* Filter police */
[name LIKE '%police%'] [@scale < 250000]{
   fill: #CDCDFE;
   stroke: #CCCCCC;
   stroke-width: 0.5;
}

/* Filter hospital/clinic */
[name LIKE '%hospital%' or name LIKE '%clinic%'] [@scale < 250000]{
   fill: #FFCCCC;
   stroke: #DCA3A3;
   stroke-width: 0.5;
}

/* Filter Educational Institutes */
[name LIKE '%school%' or name LIKE '%college%' or name LIKE '%vidya%' or name LIKE '%university%'] [@scale < 250000]{
   fill: #FFD9B3;
   stroke: #c6b29F;
   stroke-width: 0.5;
}

/* Filter Stadium/Ground */
[name LIKE '%stadium%' or name LIKE '%ground%'] [@scale < 250000]{
   fill: #DFDF9F;
   stroke: #D4D4AA;
   stroke-width: 0.5;
}

/* Name Display */
[@scale < 35000]{
   label:[name];
   font-family: Arial;
   font-size: 12;
   font-weight: bold;
   font-fill: #9C9A9A;
   font-opacity: 0.8;

   label-anchor: 0.5 0.5;

   halo-radius: 1.5;
   halo-color: white;
   halo-opacity: 0.8;

   -gt-label-auto-wrap: 80;
}

Is there any code missing? and How can I solve this Issue?


texts: Is there a tool in arctoolbox for IDistanceOp.CostPathasPolyline?

texts:
Almost all the methods in IDistanceOp have a corresponding tool in the Distance Toolset Spatial Analyst Toolbox.  However, I can't find anything in the toolbox for IDistanceOp.CostPathasPolyline.

Does anyone know if such a tool exists?


texts: Why do I have to reload the page, to show the Map correctly?

texts:
I added a map, in the way I'm supposed to, but why do I have to reload/refresh the page, to get the map (almost) as I want to? This is, what I mean:
Map before reloading the page:

Map after reloading the page:

Here is my Code:

var map;
    require([
        "esri/map",
        "dojo/dom",
        "esri/layers/ArcGISTiledMapServiceLayer",
        "dojo/domReady!"],
        function (Map, dom, Tiled) {
                ....
            var map = new Map("map", {
                logo: false
        });
});
    ...
    var luftbild = new Tiled(URL);
    map.addLayer(luftbild);

And this is how I display the map:

            ...

Karte

I'm developing with jquery mobile and the ESRI API for JS (in Visual Studio 2013)

EDIT

I managed that the map gets displayed well the first time the page has loaded:

I just wrapped the code from above with:

$(document).ready(function(/*JS code from above*/));

But I still have to reload the page, to display the map as I want to:

How can I force a reload just on this site, ASAP the site is done with loading?


texts: Difference between Layer File, Layer Package, and Map Package?

texts: What would I use them for? I know map package comprises all the layers but how is a map package different from the others?


texts: Tessellate polygon geometry in a Python script

texts:
Is there a way to tessellate a polygon geometry object, by using an existing method/function? A function call similar to feature.geometry().tessellate(). I know that there's an existing script in the toolbox which does this but I'm trying to avoid creating unnecessary layers by using processing.run("3d:tessellate,...).

And a related side question. Where can I find the existing 3d:tessellate script in my computer? I couldn't find it in C:\OSGeo4W64\apps\qgis\python\plugins\processing\algs\qgis  which is where most of the other scripts are located. I was hoping to look into the code to see how it was implemented.


texts: Finding x% of closest features to point using Network Analysis of QGIS

texts:
I have address points (Z locations), route lines, and polygon shapefiles. I need to find the route lines needed to reach 40% of all Z locations per polygon area. I have created centroids for each polygon to use as an A location. Ideally, I would want the A location to be the densest area of Z locations as I want to reach these addresses with the least amount of route miles. I then would need to do the same thing to find route lines needed to reach the closest 60%, 80%, 95%, and finally, 100% of all Z locations. All datasets have been joined with a common attribute.
I need an automated way to do this as I have 55 polygons where I need to find each 5 percentiles. I have QGIS 3.10.8


texts: Why are the organization's Maps and Layers not visible in My Content?

texts:
I'm working with a client who has an established ArcGIS Online account. They added me as an Administrator.

From the map, the Add Data button gives me access to the layers which have been configured for the organisation. But the My Content page doesn't list anything (no maps, no layers, no service definitions, etc).

Since I'm an Administrator and since I can see the layers via Add Data, shouldn't I have access to all of these items via My Content?

I can get to an individual layer's details from the map via right-click > Show Item Details. Doing so confirms that this layer is shared with the whole organisation.

Even when I choose Allow others to export to different formats on a layer, I can't see this layer when I use the Add Data From ArcGIS Online option in ArcMap. (This method works for other sites that I'm a member of.)

How can I quickly access all of the layers via My Content?

Edit: I obtained the login of another Administrator, who can see layers under their My Content page. The layers are listed as Shared amongst the Organization - but I still can't see them while logged in as myself.

There are also lots of Groups, so I added myself to all of them, using the other Administrator's login. Even so, I still can't see any items under My Content.

Are there any other settings which control who can access a layer via their My Contents page?


texts: Lat, long values extracted in QGIS do not match/overlay with polygons on map

texts:
I have been using QGIS to extract lat and long vertices from a number of shapefiles.  It worked great but for one file. When I create the lat and long for this file, the values are greater than 1,000,000.  The real lats and longs should be in the 40 to 90 and -40 to -90 range respectively.  I have tried many ideas posted in help groups but no luck. The shapefile is available at "Area-Based Conservation Measures".

If you can get the lats and longs in another software package and send them to me in a CSV file or Excel, it would be a major help. There are 250,955 features in the table.


texts:
I want to use the RQGIS

library(RQGIS)
# set the environment, i.e. specify all the paths necessary to run QGIS from
# within R
set_env()

Trying to find QGIS LTR in C:/OSGEO4~1
$`root`
[1] "C:/OSGeo4W64"

$qgis_prefix_path
[1] "C:/OSGeo4W64/apps/qgis"

$python_plugins
[1] "C:/OSGeo4W64/apps/qgis/python/plugins"

# under Windows set_env would be much faster if you specify the root path:
# set_env("C:/OSGeo4W~1")

find_algorithms(search_term = "([Pp]olygon)(centroid)")
Error in py_run_string_impl(code, local, convert) :
  NameError: name 'processing' is not defined
  attach RQGIS

I search this problem, and I knew it's python-QGIS problem.

So I tried below, but it shows this error

pip install qgis

**ERROR: Could not find a version that satisfies the requirement qgis (from versions: none)
ERROR: No matching distribution found for qgis**

What is the best way to solve this?


texts: ArcGIS Pro Calculate Field giving WARNING 002858: Certain rows set to NULL due to error while evaluating Python expression

texts:
I'm working with some county level data for the 2016 US Election.

I'm trying to calculate a blank Text field "Winner" based on two numeric fields "votes_dem" and "votes_gop" where it will take the greater of the two values in the numeric fields and input either "Trump" or "Clinton" into the new Winner field.

When I try to run my code block it returns with this warning:

  WARNING 002858: Certain rows set to NULL due to error while evaluating Python expression:   File "", line 2, in Reclass

I think there's a problem with my syntax, or just the way I'm trying to do it in general (still learning Python)

Wnner = Reclass(!Winner!)

def Reclass(Winner):
    if (votes_gop > votes_dem):
        return Trump
    elif (votes_dem > votes_gop):
        return Clinton


texts:
Geotools contains some packages dealing with KML: org.geotools.kml and org.geotools.kml.bindings. I assume it was designed to parse KML data.

Do you know how to use it?


texts:
I've recently been given a set of files to work on/analyze with QGIS, but they seem to be missing any .shp files.
Basically, the data is organized into a bunch of folders labelled land, water, etc. and in each of the folders are a few .dbf files and then a bunch of files with no extensions. Here's an excerpt from one of the ones with no extension that I opened in a text editor:

&
'                     &  è        à   ZJA   àqžVA   @bÂKA   `ShWA
-ÎYJW§·J          Q¿UJÈ»J            7YJšõ¹J
ŸWJ†Ú¹J            å [Jœ«¹J            ”ÔWJßt¹J
ŒÎUJØU¹J            ã…XJN¹J              )'YJ[¹J
ãUJƒ¹J              léWJqú¸J
ä”YJW

All the files with no extension look something like that; I've tried converting them into .shp files but QGIS still can't understand them. Do I need the .shp to do anything ?


texts: What does cumulative count cut actually do?

texts: I have a grayscale satellite image (16-bit image) which is dull (low brightness). When I apply cumulative count cut from 2 to 98 percent (Properties -> Min-Max Setting -> Cumulative count cut), it improves the visibility. I also tried to do histogram equalization but the result is not as good. So what does the cumulative count cut actually do to the image?


texts: Extract lat/long from GeoTIFF with Mercator CRS

texts:
I have a GeoTIFF file and with gdalinfo command I found the following info:

Size is 1200, 1400
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",42],
    PARAMETER["central_meridian",12.5],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-598884.437500000000000,650790.250000000000000)
Pixel Size = (996.479919433593750,-999.677795410156250)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2018:04:01 00:10:06
  TIFFTAG_DOCUMENTNAME=/datamet/data/tmp/export/RADAR/36/RRN/VMI/VMI.geotif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=IDL 8.3, Exelis Visual Information Solutions, Inc.
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -598884.438,  650790.250) (  4d32'13.20"E, 47d34'44.14"N)
Lower Left  ( -598884.438, -748758.664) (  5d56'18.16"E, 35d 4'36.70"N)
Upper Right (  596891.466,  650790.250) ( 20d26'12.34"E, 47d34'50.73"N)
Lower Right (  596891.466, -748758.664) ( 19d 2'23.69"E, 35d 4'40.94"N)
Center      (    -996.486,  -48984.207) ( 12d29'17.00"E, 41d33'32.31"N)
Band 1 Block=1200x1 Type=Float32, ColorInterp=Gray

Can you tell me a possible algorithm to associate a pair of lat/long coordinates to each point of the grid?


texts: Using ArcGIS ModelBuilder to Spatially Join using subset of point shapefile?

texts:
I would like to build a model.
The objective of the model is to subset data (based on a specific field value) from a file with 800,000 of points and spatially join (intersect) this subset data to polygon files that correspond to the time frame represented by the subset of point data.
Each polygon file represent concentrations throughout the country for a given time. So the point file is being subset for a given time frame and I need to spatially join that selection to the polygon file that corresponds to that time frame. There are 2160 polygon files.

I understand how to get the model to iterate through the point data to create the subsets. The problem I'm having is how to get the model to select the join feature for the "spatial join" tool that is appropriate for the time period represented by the point shapefile. Is there a way to use the "Get field value" tool to extract the field value from the point file as a string that can be used as a criteria to select the polygon shapefile from its folder?

I've attached a visual of the model being built thus far.
While I can do the spatial join manually, repeating the process for each hour of each day for 3 months is inefficient and introduces the potential for error.


texts: Figuring out coordinate system used in unclassified CIA document from 1983

texts:
Following unclassified document (dated 19 September, 1983) points to a non-operational historic space launch facility at these coordinates (13°40'25.0"N, 80°13'35.0"E).
https://www.cia.gov/library/readingroom/document/cia-rdp91t00712r000100290008-4
But looking above coordinates in Google Maps points to bare ground, and from public information available on this, I know the relevant facility could be either one of these two locations about 890 meters southwards.
Location one
Location two
So I wanted to know if there is way to determine geodetic system used in document. I have tried using following tool to brute force it but nothing seems to come close to real location. For example this link treats historic coordinates as EPSG:4668 (European Datum 1979) and transforms them to EPSG:3857 (WGS 84 / Pseudo-Mercator)
https://epsg.io/transform#s_srs=4668&t_srs=3857&x=80.2263889&y=13.6736111


texts: Can I reset "latlng" coordinates on Google Maps?

texts:
I have made a custom map type using the Google Maps API.

I have used maptiler to create the different zoomlayers and the 0 layer (really zoomed out) is a 256 x 256 tile. The image only covers 134 x 112 however.

I want to be able to make the map coordinates relative to the whole map and not the default earth coordinates of google maps. For example, the top-left corner should be (0,0) while the bottom-right should be (50,50).

I have received help from Heitor Chang here: https://stackoverflow.com/questions/11192281/can-i-reset-latlng-coordinates-for-a-custom-map-type

The problem lies here:

dayzProjection.prototype.fromLatLngToPoint = function(latlng) {
// pixel size of tile at zoom 0
var max = 134;

// define bottom-right corner as (50, 50)
var x = max * latlng.lng() / 50;
var y = max * latlng.lat() / 50;
return new google.maps.Point(x, y);
};

dayzProjection.prototype.fromPointToLatLng = function(pixel) {
// inverse conversion

var lng = pixel.x / max * 50;
var lat = pixel.y / max * 50;
return new google.maps.LatLng(lat, lng);
};

Can anyone help me transform these coordinates.

This may help: https://developers.google.com/maps/documentation/javascript/maptypes#Projections


texts: Order of input when calculating the difference between two DEMs over time?

texts: I am attempting to analyze coastal landscape change. I have two DEMs from significant post-storm periods for my study area. I would like to produce a simple DEM of difference using the Minus tool. My DEMs are from 2016 and 2014, if I want to see how much change occurred from 2014 to 2016 would I use 2016 DEM for Input 1 or 2?


texts:
I'm very new to GIS.

How do I find the following info?

x resolution    <--the number of real world units per pixel in the X direction
amount of translation  <-- usually 0 or a very small number
amount of rotation   <-- usually 0 or a very small number
negative of the y resolution  <-- the negative of the number of real world units per pixel in the Y direction
xMin coordinate (upper left)
yMax coordinate (upper left)


texts:
I have exported several dozen shapefiles (CRS ESPG:4326/WGS-84) from Manifold GIS  and validated the geography in QGIS. Am now editing the attribute tables in QGIS to have just two fields, an ID column and a PARCEL column.  These are mostly all single line or polygon shapefiles.

Some exports lacked a .dbf file and thus have no attribute table. Finding no method to create an attribute table in QGIS, I tested coping/pasting/renaming a .dbf file from a closely related shape export into three shapes lacking .dbf files. Then opened and edited the resulting attribute tables. The results are geographically correct with the desired attribute tables. Is this an appropriate method to create a missing .dbf file? Are there problems?

Past posts have recommended Open Office software for fixing shape .dbf files. But recent posts indicate that Open Office updates have removed that capability.


texts: Convert NetCDF(.nc) (TROPOMI Sentinel-5p) to GeoTIFF (GDAL, GRASS)

texts:
I tried converting NetCDF (Sentinel-5p) file to GeoTIFF file but I'm getting:

D:\>gdalinfo 20181119T121048.nc
Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Warning 1: 20181119T121048.nc is a netCDF file, but without any dimensions >= 2.


texts: Cross-validating kriging results using gstat R from hyperparameter

texts: I have a large dataset that doesn't require kriging (it would take too much time to compute anyways). But I'm interested to find the optimal number of samples that would be sufficient in order to apply kriging with a reasonable accuracy. What I would like to do is to compare the kriging results on the sampled datasets (N=500,750,1000,etc..) with the observed points and display a curve of the RMSE results from cross-validation. Is there a way to achieve this using gstat in R ? In scikit-learn, I would tune hyper-parameters with a grid search. N (number of samples) would be this hyper-parameter, but I don't know if I can achieve this using R. I noticed that a krige.cv function exists in gstat, would it be the way to go ?


texts:
I am using Win7.  I had QGIS 3.4.12.1 standalone installed with QGIS-OSGeo4W-3.4.12-1-Setup-x86_64.exe.

I downloaded new LTR standalone version QGIS-OSGeo4W-3.4.14-1-Setup-x86_64.exe.
This setup asks for uninstalling my previous one.
I acepted, and anfter removing previous version, new one is not present in start menu.  I can find QGIS 3.4 entry, but is empty.

I uninstalled both versions and cleaned registry form older QGIS versions but new installation have the same problem.

How should I install 3.4.14-1 version in Windows 7?

Clarifications:

Installation is ready. No error messages. Only start menu is empty.

Qgis runs when I rename qgis-ltr.bat.tmpl to qgis-ltr.bat.

But it seems to be that something is wrong with start menu. Qgis 3.4 menu is created empty and no icons in desktop. QGIS 3.4 is present in Start menu and as a directory in desktop. But both are empty.


texts: Making Geoprocessing Service results table sortable using ArcGIS API for JavaScript?

texts:
Im aware a sort() function exists, just not sure where it should be placed in my code. I have a function that draws routes on the map and i'm getting my results that are in my table from those routes. Do I put sort() after the lines are drawn, before the table is populated? or after the table is populated. Here is my code that populates my table.

myMap.Search = {
    gpGetStats: new esri.tasks.Geoprocessor(gpURL),
    getRouteStats: function () {
        console.log(dojo.byId('resultsLayer').value)
        var params = {
            "ZIP_CRID": dojo.byId('resultsLayer').value
        };
        resultsLayer.Search.gpGetStats.execute(params, resultsLayer.Search.handleStats);
    },

    handleStats: function (results, messages) {
        console.log(results);
        // content.sort();

        var content = "";
        content += '
';
        content += '
 ZIP CRID 	 Bus 	 Res 	 Tot 	 &lt;200 	 DS Key
';

        dojo.forEach(results[0].value.features, function (crid, i) {
            console.log(crid.attributes.ZIP_CRID + ', ');
            content += '
' + crid.attributes.ZIP_CRID;
            content += '	 ' + crid.attributes.BUS_CNT;
            content += '	 ' + crid.attributes.RES_CNT;
            content += '	 ' + crid.attributes.TOT_CNT;
            content += '	 ' + crid.attributes.LT_200_IND;
            content += '	 ' + crid.attributes.DS_KEY + '
';
        });
        content += '
';

        dojo.byId('details').innerHTML = content;
    },
}

Also a Dojo Grid would be acceptable


texts: Plotting Geopandas Gallery example fails?

texts:
I have been having trouble getting the figures from the example provided in Geopandas Gallery examples. When I ran these scripts in Jupyter console:

import geopandas
import geoplot
import matplotlib.pyplot as plt

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
africa = world.query('continent == "Africa"')
ax = geoplot.cartogram(africa, scale='pop_est', limits=(0.2, 1), edgecolor='None', figsize=(7, 8))
geoplot.polyplot(africa, edgecolor='gray', ax=ax)

I will only get the this output:

<matplotlib.axes._subplots.AxesSubplot at 0x7f31fff31d30>

I have tried including the %matplotlib inline at the beginning of my Jupyter console, and I also have tried plt.show() at the end of the script. However, my console did not display the figure. I am looking into getting this figure:

Does someone have an idea to visualize the output a matplotlib.axes._subplots.AxesSubplot object?

I am on an Anaconda bundle using:
Python 3.6,
Geopandas 0.5.1,
Matplotlib 3.1.1


texts: How can I join several lines after offset in Python

texts:
I had a layer, and I did one offset of it. But the problem is that now I have a offset but the lines are separated.

How can I join them in Python? Because I read about several function as snap, snapped line, snapped geometry, merge, etc.

I wanted to know how can join theses lines, but with Python code.

In the above image you can see the lines after offset, and the next was the code.

import processing

layer = iface.activeLayer()

parameters = { 'DISTANCE' : 5, #offset of 5 meters
               'INPUT' : layer,
               'JOIN_STYLE' : 1,
               'MITER_LIMIT' : 1,
               'OUTPUT' : 'TEMPORARY_OUTPUT',
               'SEGMENTS' : 8 }

processing.runAndLoadResults('qgis:offsetline',
               parameters)


texts: (pyqgis) How to get temp object created by the processing tool "qgis:saveselectedfeatures"

texts:
Usually when I try to get object:

buffering = processing.runalg('qgis:fixeddistancebuffer', in_layer, 100, 10, False, None)
buf_layer = processing.getObject(buffering['OUTPUT'])
QgsMapLayerRegistry.instance().addMapLayer(buf_layer)

But when I do the following it does not work:

saving = processing.runalg('qgis:saveselectedfeatures', in_layer, None)
output_layer = processing.getObject(saving['OUTPUT'])
QgsMapLayerRegistry.instance().addMapLayer(output_layer)

The error is

KeyError: 'OUTPUT'

Anyone knows how to get this object?


texts: Use a Tiledmap service as basemap instead of an Arcgis.com map in AMD style coding

texts:
How will I define a arcgisTiled map service as the base map (instead of “topo”) as against using arcgis.com map services while defining a new map in AMD style?

map = new Map("map", {
          basemap: "topo",
          center: [-122.445, 37.752],
          zoom: 14
        });

I would like to replace "topo" with a map tiledservice from my server:

http://ourserver/arcgis/rest/services/tiledmapservice/MapServer

Any thoughts?


texts: Arcpy python field calc syntax error variable creation

texts:
I am fairly new to python. I'm trying to set this function into a variable, but i can't seem to get it to work.

Can you take a look at my script?

import arcpy
arcpy.env.workspace = "M:\Parks\HQ\OPERATIONS\CPI\Asset_Management\GNSS Projects\Parks Inventory\MergedParks\Ram_Falls_PP\CW_Ram_Falls_PP\CW_Ram_Fall_PP.gdb"

a2 = """def rep_field(in_fld, rep_value):
            targets = ['W:\OPERATIONS\CPI\Asset_Management\GNSS Projects\Parks Inventory\Ram Falls PP\RAM FALLS PP JULY 28 2015 JR~files','W:\OPERATIONS\CPI\Asset_Management\GNSS Projects\Parks Inventory\Ram Falls PP\RAMFALLSPP_JULY28_2015_JL~files']
            for targ in targets:
             in_fld = in_fld.replace(targ, rep_value)
            return in_fld"""

a1 = """rep_field (!Picture!,'\\\env.gov.ab.ca\Parks\CityWorks\Referenced\Ram_Falls_PP')"""

a11 = """rep_field (!Picture1!,'\\\env.gov.ab.ca\Parks\CityWorks\Referenced\Ram_Falls_PP')"""

arcpy.CalculateField_management(in_table="Bench_and_Table", field="Picture", expression = a1, expression_type="PYTHON_9.3", code_block = a2)
arcpy.CalculateField_management(in_table="Bridge", field="Picture1", expression = a11, expression_type="PYTHON_9.3", code_block = a2)

It is now giving me a

  ERROR 000989: Python syntax error: Parsing error IndentationError: unindent does not match any outer indentation level (line 3)
  Failed to execute (CalculateField).


texts: QGIS Geometry Generator - Connect Point Data with Lines

texts:
I have two point data sets; 'Person' (in red) and 'Connections' (in blue).
They have matching columns 'Person_UID'.
I'm trying to connect the red points to all the blue points that have matching Person_UIDs.
So far I am only able to make one connection for each UID in'Connections' the point that was first digitized.

Code:
make_line(
  $geometry,
  geometry(get_feature('Connections', 'Person_UID', "Person_UID"))
)

I'm using QGIS 3.12.0

Screen grab attached.


texts:
For example, I have 2 points and I want to know if they intersect. If yes, I need to delete first points. Should I use OpenLayers .intersects method? Because it greatly depends on zoom level.

Or should I use distanceTo method?

Or how should I do control intersections?


texts: Automating specific workflow using QGIS?

texts:
For a project of mine, I'd like to solve a spatial analysis using QGIS 2.8.

Until now, I have never automated anything GIS-related except for some rather easy steps in ArcGIS using the ModelBuilder.

This is what's supposed to be done:

Select all polygons from shapefile A that overlap a similar shapefile B (-> Select by location)
Extract selected polygons from A creating shapefile C, B serves as clip feature (-> Clip)
Add a columns reading area for each of the extracted features of C (-> Add geometry columns)
Calculating the percentage of each feature's (C) area in regard to all extracted polygons (= sum), store this information in a new column (I'm actually not sure how to do this except manually using the individual .dbf-file, the GroupStats plugin works fine but cannot write  the newly calculated information in the attribute table)
Calculate a new value in a new field in C's attribute table based on the area's percentage (see 4) and an existing value (derived from A, see 2) (-> Field Calculator)
Summing up the calculate values (see 5) in a new field (-> Field Calculator)
Merge all selected polygons from C (-> Editor or Dissolve) (see 2)

For me, the major problem is that this workflow needs to be done some hundred times at least. Therefore, if this can be easily done using the Modeler, the Python console or something (my Python knowledge is somewhat limited but not entirely non-existent), I am willing to dedicate some time to delve into tutorials or instruction manuals in order to be better off next time.

Here's my latest model. How can I perform the fourth task?

Trying to add underdark's script from file resulted in this configuration error message:

  Traceback (most recent call last):
    File "/home/tjk/.qgis2/python/plugins/processing/gui/ProcessingToolbox.py", line 221, in executeAlgorithm
      action.execute()
    File "/home/tjk/.qgis2/python/plugins/processing/script/AddScriptFromFileAction.py", line 45, in execute
      settings = QSettings(QSettings.NativeFormat, QSettings.UserScope, 'QuantumGIS', 'QGis')
  NameError: global name 'QSettings' is not defined

I tried the solution to https://stackoverflow.com/questions/18769383/qgis-preferences-in-a-python-plugin but nothing happened.

Shouldn't QSettings be defined automatically while I'm using the software, saving my personal preferences?


texts:
I am being baffled by how apparently poorly NetworkX reads a shapefile and builds a graph out of it.

Below is a graphical example of a fake network built with 27 polylines all snapped together, so there are no topological errors (ArcGIS 10.3.1, License: Advanced). The coordinate system is the British National Grid (OSGB36).

Now, when doing

import networkx as nx
import matplotlib.pyplot as plt

G=nx.read_shp('C:\Users\MyName\MyFolder\TEST.shp')
nx.draw(G, with_labels=False)

The result is a graph plotted according to the spring layout, so the shape of the file is not preserved. Also, the image shows three disconnected components, which are absolutely fabricated, since all segments in my shapefile were snapped. Also, the coordinates of the nodes are imported as labels.

I wonder what the usefulness of the nx.read_shp() function is, if you are not allowed to reproduce the shapefile exactly. This is especially true if you want to make sure that the connections are recreated accurately.

This leads me to my question: how to play with NetworkX so that the graph is not plotted in some meaningless way, but reproduces the network shown in the shapefile?


texts: Any way to style LineString-features with Stroke and Fill?

texts:
I am using the latest version of OpenLayers to draw some features within a map rendering data from GeoJSON-files. Within each of these files I have a FeatureCollection containing a couple of LineString-features.

So far the following style is applied, which works fine and draws each feature as a thin dark line:

return new Style({
  stroke: new Stroke({
    color: '#333333',
    width: 1
  })
});

I know, that I can not use Fill to style LineString-features (if I am right). However, I really want to style these lines with two colors using a inner line (1px;red) and a border/stroke around this inner line (1px;blue).

Is there any possibility to accomplish this?

A first solution would be to copy the features drawing the lower ones blue and 3px wide and the upper ones red and 1px wide - but this would also double the amount of features.


texts: Creating fishnet from thematic polygon data?

texts:
I need to create a regular grid containing attributes from other thematic data from ecological data. I am in effect creating a homogenous polygon raster and adding information from another dataset which contains irregular-shaped polygons to the regular shaped raster cells . I need to therefore clip or union the fishnet with the thematic data - but the only problem is that inside the fishnet cell only one value is permitted (because this will eventually become a raster-formated dataset most likely a geotif. How do I decide which value the cell obtains? - some features are more important than others and should be weighted, other cells would just obtain the value from the polygon with the largest area inside the cell or possibly because it is location centrally inside the cell.

Below is an illustration of the data in ArcGIS.


texts: Joining postcodes to coordinates in a .csv file

texts:
I have a .csv of land registry sales data that includes full addresses and postcodes.

I'm trying to geocode all the data so that I can plot it in QGIS but my issue is that I have 2 million records that I'm not sure I can feasibly run through any commercial geocoders without paying heaps (if anybody has any ideas let me know)

Instead I've downloaded Ordinance Survey's code point dataset which lists coordinates for all of the postcodes in Britain. I'm hoping that there is some way that I can join those records with those in my landreg dataset so that I can have some degree of accuracy in plotting my points.

Any idea on how to do it?


texts: Updating/refreshing/redrawing vector layer in OpenLayers?

texts:
I am creating a layer groups with:

      var vectorSource = new ol.source.Vector({
        features: (new ol.format.GeoJSON()).readFeatures(geoJson, {dataProjection: 'EPSG:4326', featureProjection: 'EPSG:3857'}),
      });

      var vectorLayer = new ol.layer.Vector({
        source: vectorSource,
        style: observationStyle
      });

      vectorLayer.set('name', observationGeoJson.features[i].properties.report_id);

      layerGroup.push(vectorLayer);

      ...

      var newLayerGroup = new ol.layer.Group({
        layers: layerGroup,
        style: style
      });

The GeoJson data source is then modified (the position of a point is changed). What is the best way to update the layer group? I have tried map.updateSize();, removing and adding the layer, .getSource().change() and .getSource().clear() + getSource().addFeatures(newFeatures)
But the map is always rendering both the old geometry and the new geometry.


texts: Most 3D visualisation solutions I've seen are drapes of imagery over a surface like a TIN.  Is there any open source software out there that can permit viewing of 3D data (e.g. Leica ADS data) and creation of a terrain following vector layer similar to what the Stereoanalyst plug-in for ArcGIS does and which doesn't cost a lot?


texts:
I am a software engineer that is very new in the geo world. So excuse my ignorance if this question seems simple or if terms I use are totally wrong.

Problem:

I have been tasked with taking multiple geotiffs and serving them from a map tile system on our servers. From my research so far, I believe that I need to convert these geotiffs to mbtiles. I've got the serving of the mbtiles part figured out, what I really need help on is coming up with a quick way to convert the geotiffs into mbtiles.

Things I've tried:

I currently have a script that can do this conversion of 56GB of Geotiffs in about 7.25 hours using gdalbuildvrt, gdalwarp, and gdal2tiles. Here is what it does:

ls -1 *.tiff > files.txt
gdalbuildvrt -input_file_list files.txt project.vrt
gdalwarp -t_srs EPSG:3857 -of VRT project.vrt project_mercator.vrt
gdal2tiles.py project_mercator.vrt project.mbtiles

It works but I'd like to improve on this.

I've heard that it might be faster to use TileMill to export, but so far my attempts to import the geotiffs into TileMill have failed. I have 158 geotiffs and I don't want to have to "Add Layer" for each one using the interface.

I thought I might be able to import the .vrt file I generated in my own script but when I try to add it in TileMill I get an error that says, "GDALWarpOptions.Validate() hSrsDS is not set".

Question to the community:

What is the best way to import multiple geotiffs into a TileMill project?


texts:
I want this layer...

to have the same division/cut/split like this..

I tried the Split Features tool but this software (qgis 2.16.0) would hang then eventually force close.

Any method/tool that you can suggest?


texts: Workstations versus desktop computers for GIS applications

texts:
ArcGIS for Desktop has native 32-bit data processing limitations, which is somewhat offset by 64-bit background processing capabilities.  Additionally, multiprocessing and parallel processing capabilities are limited on ArcGIS for desktop.

Given these software limitations, what are the pros and cons of using a workstation over a desktop computer for geoprocessing with ArcGIS for Desktop?  For example, how do vector and raster based processing speed tests compare between these classes of computers?  Is the additional expense of a workstation worth the cost of the software bottleneck?


texts: what are the main advantages of using a Double Stereo projection?

texts: I'm not per say a GIS guru, I'm more of a user of the technology. On a recent project we were asked to analyze, provide, and submit all our analysis using a Double Stereo projection. What are the main advantages of using it vs. UTM, or any other projection.


texts: Getting all (or selected) values from single field in table using ArcPy?

texts:
I'm looking for a way to get all of the values from a single field in an ArcGIS table.  Reasons for doing this might include pasting the values into a text document or SQL where clause.

Since ArcGIS 10.1 does not seem to have a way to do this (like right-click the field > Copy), I would like to see if I could write an ArcPy script to do this for me.  The following code will write the values to the Python window if I pass it an explicit tablename and fieldname:

import arcpy
cur = arcpy.SearchCursor("myTableName")
fieldname = "myFieldName"

for row in cur:
    print row.getValue(fieldname)

del row, cur

However, to make this a dynamic tool, I would like to determine a way to identify:

the active table (in the table window), and
the active (highlighted) column/field in that table

at which point I can dynamically change the tablename and fieldname in the code above, based on what is selected.

Any ideas?


texts:
I'm trying to split a shapefile of geometry type 'line' at a shapefile of geometry type 'point. These points are snapped to the lines and are the start points/end points of other line shapefiles.

Something like this :
The points are snapped to the black lines.

As of now I'm working with Geopandas/Shapely/Fiona.
I've searched for a native function in Shapely that performs the split; however, I've had no luck in finding a solution. This is also true when searching through Shapely's official docs.
This is odd since the split is among the basic geoprocessing tools.

My goal is to give a function of the two shapefiles (lines & points) and split the lines at points (like ArcGIS' tool). The split lines should then be stored in a different shapefile.
This seems easy to do but the lack of documentation is a bump in the road.

Any ideas on how I can do this in python using Geopandas/Shapely?


texts: How to make icon marker size variable depending on data value using CartoDB?

texts:
Is there any way in CartoDB where I can make the size of individual icon markers, visualized on a CartoDB map, variable and set according to its data value in a given data range?

Similar to the CartoDB "Bubble" visualization but i would like to be doing this for other Icon markers other than circle markers.

Anyone having any idea I am most grateful for any clue or maybe recommendation for other tool than CartoDB that can do this if CartoDB cannot do it.


texts: Updating Z values on points using UpdateCursor

texts:
I have created the following script, which moves the points from my_input_points_fc to my_output_points_fc by user input
Furthermore, the z-values in my_output_points_fc should also be updated. However, after executing the script I get zero values in the my_output_pointy_fc under the height_field.
What am I doing wrong?
import arcpy
arcpy.env.overwriteOutput = 1
arcpy.env.workspace = "my_path"
inputFC = "my_input_points_fc"
outputFC = "my_output_points_fc"
xValue = -1
yValue = -1
while xValue == -1:
    try:
        xOffset = input("Please enter the X coordinate by which the point should be move: ")
        xValue = float(xOffset)
    except Exception:
        print()
        print("The entry was not correct! Please try again")
while yValue == -1:
    try:
        yOffset = input("Please enter the Y coordinate by which the point should move ")
        yValue = float(yOffset)
    except Exception:
        print()
        print("The entry was not correct! Please try again")
spatialRef = arcpy.Describe(inputFC).spatialReference
arcpy.management.CreateFeatureclass(arcpy.env.workspace, outputFC, "POINT",inputFC, has_z = "ENABLED", spatial_reference= spatialRef)
cur = arcpy.da.InsertCursor(outputFC, ("SHAPE@XYZ"))
cur2 = arcpy.da.SearchCursor(inputFC, ("SHAPE@X", "SHAPE@Y", "height_field"))
for row in cur2:
    insert_row = (row[0] + xValue, row[1] + yValue, row[2])
    cur.insertRow([insert_row])
print ("Points have been moved")
del cur
del cur2
del row


texts: Python Scripting in QGIS: attribute index out of bounds [0;0]

texts:
I am currently writing an Python Script within QGIS to connect points of two specific point layers:

Attributes of point layer 1:

Attributes of point layer 2:

For each point in point-layer1, the python script loops over all points of point-layer2 and compares the attribut "direction" (of point-layer1) with the attribut "location" (of point-layer2). If they are equal, a new line should be created and stored into the output line-layer (together with the lineID (a line counter) and the direction name + location name).

Here is my Python Script:

# name : ConnectPointsByValue_script
# date : 27.12.2017
# description: connects points from evis with points of location_coordinates
#       if the attributes "direction" and "location" are equal

from PyQt4.QtCore import QVariant
from qgis.core import *
from PyQt4.QtCore import *

#***********Constants
#
LINEID = 'lineID'
POINTID1 = 'pointID_Location'
POINTID2 = 'pointID_Direction'
#
#***********

# input parameters
##Counters=vector
##Locations=vector

inputLayerCounters = processing.getObject(Counters)
inputLayerLocations = processing.getObject(Locations)

crs = inputLayerCounters.crs().toWkt()

# Create the output layer
outLayer = QgsVectorLayer('Linestring?crs='+ crs, 'connector_lines2' ,
'memory')
outdp = outLayer.dataProvider()

#get the field index
outLineID = outdp.fieldNameIndex(LINEID)
outID1 = outdp.fieldNameIndex(POINTID1)
outID2 = outdp.fieldNameIndex(POINTID2)
lineIDcnt = 0

# loop through points of iCounters and connect each point with the matching
location point
for feat1 in inputLayerCounters.getFeatures():
   geom1 = feat1.geometry()

   #loop all points from iLocations
   for feat2 in inputLayerLocations.getFeatures():
      geom2 = feat2.geometry()

      # get the feature's attributes
      attrs1 = feat1.attributes()
      idx1 = inputLayerCounters.fieldNameIndex('direction')
      attrs2 = feat2.attributes()
      idx2 = inputLayerLocations.fieldNameIndex('location')

      #check if attributes are equal
      if feat1.attributes()[idx1] == feat2.attributes()[idx2]:
          lineIDcnt +=  1

          #if true: create line and add to output line layer
          connector = QgsFeature(outLayer.pendingFields())
          connector.setGeometry(QgsGeometry.fromPolyline([geom1.asPoint(),
                                                    geom2.asPoint()]))

          #record the IDs
          res = connector.setAttribute(outLineID, lineIDcnt)
          res = connector.setAttribute(outID1, feat1.attributes()[idx1])
          res = connector.setAttribute(outID2, feat2.attributes()[idx2])

          #add the feature to the layer
          res = outdp.addFeatures([connector])

 # Add the layer to the Layers panel
 QgsMapLayerRegistry.instance().addMapLayer(outLayer)

 print "Done .. "

If I execute this script, I get an error which says:

attribute index -1 out of bounds [0;0]

I think I am doing something wrong with indexing the attributes of my points. Any suggestions on how to solve this?


texts:
I have created this array:

array = ([arcpy.Point(-103.85769159, 32.181106371), arcpy.Point(-103.85769159, 32.182755629), arcpy.Point(-103.85575241, 32.182755629), arcpy.Point(-103.85575241, 32.181106371)])

and would like to create a polygon feature from it:

polygon = arcpy.Polygon(array)

My map units are in Decimal Degrees and I'm using NAD27, no projection.

Whenever I try the line to create the polygon, I get an error:

"CreateObject cannot create geometry from inputs"

Why is this happening?  It seems so straight forward.


texts: pyQGIS returns invalid layer with QgsVectorLayer in QGIS 3

texts:
I have the following problem :
I've installed all pyQGIS environment on Windows 7. I work with 64 bits version of Python 3.6 and QGIS 3.0 and PyQt5. My Windows paths are ok.

My objective is to create an empty QGS project with stand-alone python script and add layers into.

But when I try to read a shapefile (the shp is good when I open it in QGIS GUI) I get an invalid QgsVectorLayer :

import os

from qgis.core import QgsApplication, QgsProject, QgsVectorLayer
from PyQt5.QtCore import QFileInfo

QgsApplication.setPrefixPath("C:\\QGIS 3.0\\apps\\qgis", True)
qgs = QgsApplication([], False)
qgs.initQgis()

pathQGS = "C:\\path\\to\\project.qgs"
layerPath = "C:\\path\\to\\layer.shp"
fileInfo = QFileInfo(layerPath)
path = fileInfo.filePath()
baseName = fileInfo.baseName()

project = QgsProject.instance()
if os.path.isfile(pathQGS):
    os.remove(pathQGS)
project.setFileName(pathQGS)
layer = QgsVectorLayer(path, baseName, 'ogr')
if not layer.isValid():
    print("Layer is not valid")
project.addMapLayer(layer)
project.write()

qgs.exitQgis()

How can I solve this ? I checked many other topics without finding any answer to my problem.

Can you please help me ?


texts:
Thanks to Best Landsat-5 TM band combination for detecting fire scars I am able to calculate burnt area in an absolute way (Normalized Burn Ratio).

Does a ratio or delta operation exist to calculate the regrowth area between two scenes?


texts: Constructing least-cost path between many villages and nearest available road using ArcMap and Network Analyst

texts:
I am using ArcMap.
I am interested to construct a least-cost path between villages and the nearest available road by using land gradient and location of rivers to create a cost raster. Hence,  I am interested to construct the least-cost paths between many points (the number of villages is 1584) and a line (nearest road). However, most of the tutorials are for only two points (one starting and one destination point).
Any suggestions?


texts: Problem with pushButton handling - pyqgis

texts:
I have made a pushButton - new window relation. Now, when I'm trying to 'link' methods with pushButtons in new window they don't react.

Here is what I wrote in metod, which opens new window in plugin:

def handleButton(self):
    self.sc = PluginConfigDialog()
    self.sc.show()
    result = self.sc.exec_()

I have also added this line in __init__:

self.dlg.pushButton.clicked.connect(self.handleButton)

But when I'm trying to add a alike action with a pushButton_2 in new window, to open another one (third window) nothing happens. Here is code I wrote:

def handleEdit(self):
    self.s = PluginEditDialog()
    self.s.show()
    result = self.s.exec_()

I also added this line to __init__:

self.d.pushButton_2.clicked.connect(self.handleEdit)

I have also added imports for this new windows and create the dialogs (after translation) and keep reference in __init__

Do you know what I'm doing wrong?


texts: Adding gdal nextgis borsch as a cmake external project

texts:
I want to add GDAL as a dependency in a modern cmake project using a superbuild pattern. Since finding a cmake aware project sounded easier, I have been trying to add https://github.com/nextgis-borsch/borsch as an external project, but without success.
The most likely hypothesis is that something is wrong in my call to external project. Here is the external/upstream/gdal/CMakelists.txt that I'm using:
cmake_minimum_required(VERSION 3.12)

find_package(GDAL QUIET)

if(GDAL_FOUND)
  get_property(_loc TARGET GDAL::GDAL PROPERTY LOCATION)
  message(STATUS "Found GDAL: ${_loc} (version ${GDAL_VERSION})")
  add_library(gdal_external INTERFACE)  # dummy
else()
  message(STATUS "GDAL could not be located, downloading and building version 3.1 instead.")

  include(ExternalProject)
  ExternalProject_Add(gdal_external
      URL                   https://github.com/nextgis-borsch/lib_gdal/archive/v2.4.0.tar.gz
      UPDATE_COMMAND        ""
      CONFIGURE_COMMAND     ""
      BUILD_COMMAND         ""
      INSTALL_COMMAND       ""
      DOWNLOAD_NO_PROGRESS  1
      LOG_CONFIGURE         1
      LOG_BUILD             1
      LOG_INSTALL           1
      CMAKE_ARGS            -DCMAKE_INSTALL_PREFIX=${STAGED_INSTALL_PREFIX}
    )

    set(
      GDAL_ROOT ${STAGED_INSTALL_PREFIX}/gdal
      CACHE PATH "Path to internally built gdal installation root"
      FORCE
      )
    set(
      GDAL_INCLUDE_DIRS ${STAGED_INSTALL_PREFIX}/gdal/include
      CACHE PATH "Path to internally built gdal include directories"
      FORCE
      )
    set(
      GDAL_LIBRARY ${STAGED_INSTALL_PREFIX}/gdal/lib
      CACHE PATH "Path to internally built gdal library directories"
      FORCE
      )

endif()

I think I am incorrectly building the project, but the documentation of the nextgis-borsch project did not help me figuring out my mistakes, as it seems oriented towards an expert audience (what I'm not).


texts:
I have a shapefile of the 8830 election precincts in Texas as of the 2016 election.  I also have shapefiles of the Texas US House districts and the Texas Senate and House districts.  I want to generate tables that match the election precincts with the districts that they are contained in.

Using the functions available (Within, Contains, Touches, etc) all see to either produce too many rows, or too few.  One problem I identified is when an election precinct is matched with a district because they share a border, although the election precinct is in another district.  As far as I know, no election precincts cross any district border.

I got this far using QGIS and SpatiaLite.  I'm a newbie at GIS.  This is my first project.  Apologies if I have left out any critical information.


texts: Is there any extension to "optimize" the creation of maps in QGIS?

texts:
is there any extension to "optimize" the creation of maps in QGIS?
i searched the repository, but found nothing useful, maybe i searched in a wrong way..

there are some things i am really missing:

rulers to place the map-elements,
the possibility to use thousand-seperator,
the possibility of changing the elements propertiers by right-click-menu or at least in an easier way as at the moment would be pretty nice.


texts: Why does the ArcGIS Sample tool skip rows and points?

texts:
I'm using ArcGIS 10.4 and want to export XYZ values of a raster (type: GRID) to a text file. I'm following the procedure recommended by ESRI support. However, it seems that the Sample Tool in spatial analyst toolbox skips a bunch of points. When I display the XY-Data certain rows and points are incorrect, either 0 or -9999 (not displayed).

I tried this with different input data, and i cannot find an explanation at all. I thought this might be because of the ArcGIS instance i was runnning, but I tried the sample tool in an freshly started ArcCatalog window. Almost the same happens there, just that the points are not 100% identically.

Is this a bug, or am I doing something completely wrong? How do I fix this?

It seems to be related to this question from 3 years ago.

Edit: Several new tries (see suggestions below), resulting in different patterns. Could still not find the solution to this problem.
Another screenshot with nearly the same area, different zoom, different color scale. -9999 is bright yellow, 0 is dark color.


texts: How can a statistically significant difference between rasters be calculated?

texts:
so after some gruelling use of MaxEnt I have four occurrence probability rasters for the four quarters of the year, each one is about 70MB+. I want to test whether they are significantly different from each other. I have been told that a Mantel test is the best option, however it seems to only be able to deal with two rasters at a time, and with such large rasters and numbers of permutations R has been failing me in terms of memory:

rsession(14001,0x7fff7be68000) malloc: * mach_vm_map(size=153281630035968) failed (error code=3)
 error: can't allocate region
 set a breakpoint in malloc_error_break to debug
rsession(14001,0x7fff7be68000) malloc:  mach_vm_map(size=153281630035968) failed (error code=3)
 error: can't allocate region
* set a breakpoint in malloc_error_break to debug
Error: cannot allocate vector of size 142754.6 Gb
In addition: Warning message:
In is.euclid(m1) : Zero distance(s)

I wonder if there is some way of doing a Mantel test or similar within QGIS, or whether anyone can help?


texts: Calculating climatological albedo using MODIS satellite images?

texts:
My question is based on the paper Using NOAA AVHRR and SPOT VGT data to estimate surface
parameters: application to a mesoscale meteorological model.

The authors have calculated albedo, surface emissivity and thermal inertia for a small region in Europe using the AVHRR data.

In a similar way I have land cover information at 30 meter resolution derived from Landsat images.

I need to calculate  albedo for a small region and a climatological period(say Boreal summer). I used the Semi Automatic Classifer Plugin from QGIS to download Landsat images for the period June 1st through August 31st(2014 and 2013) and I got very few images over my region. I  need at least 50 to 70 images(more would be better).

I believe  MODIS satellite images provide an alternative to  calculate the albedo. It is at a lower resolution but I believe a downscaling can be done. My question is can it provide sufficient number of passes in order to come up with a climatological albedo ? The presence of clouds is not a problem as cloud masking can be done.

My ultimate purpose is to use these parameters such as albedo in a meteorological context such as with a weather model that predicts the weather.


texts: Solar Diffuse Radiation Model on an Inclined Surface

texts:
I have implemented the Perez 4 (as numbered by D., Yang 2016 ) diffuse radiaition transposition  for a solar PV model. It seems to work but I have no way of knowing if it really is.
My panels have a vertical tilt, are at a latitude of -79.5°. The Perez 4 model produces peaks of  1.9 and a mean value of 0.91 where an isotropic model would have a constant sky view factor of 0.5.
Is this reasonable or is there something wrong with my implementation? I had not expected the mean of the Perez to be so much higher than the isotropic.


texts: Interpolate data from multiple gpx tracks

texts:
I have a number of GPX tracks of the same road. But due to GPS error they are not accurate enough (20-30m off).

I'm wondering if I could interpolate those multiple gpx tracks into single one. If an error is random, interpolation should give me  more accurate data.

Are there any algorithms to interpolate multiple (non-ordered) points?

Update I've found the correct search term. It's map inference gps traces in google scholar . There you can find "cited by" and related articles.


texts: Formulation of CRS (EPSG2438 to EPSG4326) conversion in Excel

texts:
I have about 14000 points with coordinates EPSG2438 and want to convert the coordinates into WGS84 (EPSG4326). I want to know that is there a way to make Excel formula to convert these coordinates. In QGIS (and I suppose in other GIS packages as well) when we want to change CRS, we can see certain information associated with each EPSG. For Example for EPSG4326

Extent: -180.00, -90.00, 180.00, 90.00
Proj4: +proj=longlat +datum=WGS84 +no_defs

and for EPSG2438

Extent: 121.50, 28.22, 124.50, 53.56
Proj4: +proj=tmerc +lat_0=0 +lon_0=123 +k=1 +x_0=500000 +y_0=0 +ellps=krass
+towgs84=15.8,-154.4,-82.3,0,0,0,0 +units=m +no_defs

Can this associated information be used to convert from EPSG2438 to EPSG4326 and vice versa?


texts: Loading offline tiles as QgsRasterLayer in C++ or Python?

texts:
I have created tiles via Qtiles plugin in QGIS.

Here is the code snippet I have so far:

QString url = "file:///C:/Users/USER/Desktop/TiledMap/{z}/{x}/{y}.png";
QgsRasterLayer * lyr = new QgsRasterLayer ( url, "tiledLayer", "wms" );
QgsMapLayerRegistry::instance()->addMapLayer(lyr , true);

The layer does not show up though. What am I missing? Is there something wrong with the URL?

I originally created the tiles via QTiles plugin using the following options:

Minimum zoom level: 5, Mazimum zoom level: 12


texts: Calculating distance between matching intersections using QGIS?

texts:
I have a streetmap with the coordinates of some speed limit signs (which mark either the beginning or the ending of a speed limit and hold the applicable speed limit in km/h). With regard to the direction of driving, the signs are always positioned on the right side of the road:

What I want to do now is to create a new line layer which shows the speed sections (the length in km of the speed section and the speed limit in km/h itself).

As you can see the red rectangle in the image above, I will cover the marked area in an example:

In order to do that, I guess I will have to create normal vectors with the points of the speed signs and the streets itself in order to get intersections directly on the streetline. I know how to do that by now. But how can I "drive" along the streetline to sequentially check each speed sign in order to extract the speed sections? Let's take the picture above as an example: If I start south at the beginning of the 70 km/h speedsection, I have to drive north (because it is the driving direction of the street) and check each speed sign along the line for the end of the speed section. After that, I search for the next beginning speed section and so on.

How can that be realized?
Please comment if you need more details on that issue.


texts: Adding StreetView add into Google MyMaps

texts: A customer uses Google MyMaps to track relevant location data. They want to be able to use the StreetView addin in their map but I can't seem to find a way to implement it.


texts: Sorting by shape token in ArcPy cursors?

texts:
I have a script that uses the sql orderby parameter to sort my arcpy cursors by a user defined field or fields and deletes duplicates, but keeps the duplicate with the highest value. I'd like to be able to sort by geometry to identify/group features that are spatial duplicates and delete those that are duplicates (but keeping the duplicate with the highest value in another chosen field which is why I can't just use the Delete Identical tool).

I have tried several ways of doing this, but as I understand, shape tokens such as SHAPE@, SHAPE@WKT, etc. cannot be used in the sql clause. Is there another way to do this or a workaround? I also tried to add a field and fill with the SHAPE@WKT, but it didn't fill the field.

Here is the code I've got that works for other fields, other than the geometry field.

import arcpy
from itertools import groupby
from operator import itemgetter

input_layer = arcpy.GetParameterAsText(0)
case_fields = arcpy.GetParameterAsText(1)
max_field = arcpy.GetParameterAsText(2)
case_fields = case_fields.split(";")
case_fields = [str(x) for x in case_fields]

sql_orderby = "ORDER BY {}, {} DESC".format(",".join(case_fields), max_field)

with arcpy.da.UpdateCursor(input_layer, "*", sql_clause=(None, sql_orderby)) as cursor:
    case_func = itemgetter(*(cursor.fields.index(fld) for fld in case_fields))
    for key, group in groupby(cursor, case_func):
        next(group)
        for extra in group:
        cursor.deleteRow()


texts: Delete duplicate geometry in postgis tables

texts:
After - I don't know what happend - all my entries in my PostGIS tables are doubled! I tried this to delete them but it does not delete any/all duplicates:

DELETE FROM planet_osm_point
       WHERE osm_id NOT IN (SELECT min(osm_id)
                        FROM planet_osm_point
                        GROUP BY osm_id)

or this:

DELETE FROM planet_osm_point
WHERE osm_id NOT IN (
    select max(dup.osm_id)
    from planet_osm_point as dup
    group by way);

EDIT:

I finally found an easy way, which is working in my case:

DELETE FROM planet_osm_point WHERE ctid NOT IN
(SELECT max(ctid) FROM planet_osm_point GROUP BY osm_id);

found on this page: http://technobytz.com/most-useful-postgresql-commands.html


texts: In QGIS I have a polygon SHP layer. The layer contains rings. I am  reshaping the polygon layer by cutting away large areas of the polygons. When cutting away a part of the polygon that contains a ring, what happens to the ring? Is it completely removed?


texts: Adding multiple pages in QGIS print composer?

texts:
Trying to add a separate page in QGIS 3. Since the composer was changed in QGIS 3.x.x I can't find the field for multiple pages that is listed here How to add pages in print composer.

Where did it move?


texts: Issues translating 16-bit Int (NetCDF) to GeoTIFF with gdal_translate

texts:
I am attempting to translate the .nc dataset to GeoTIFF but run into issues with the colors of the final TIFF file. It is a 1-band dataset (lat, lon, temperature (K).

Command I am using:

gdal_translate -of GTiff LIIN50_202002060900.nc whydoesntthiswork.tiff

This runs without complaint, but the TIFF produced appears over saturated. ie. the values where temperature is higher are far to bright. This is compared to a plot of the NetCDF in Python.


texts: Setting Spatial Reference for Geometry Bag using ArcObjects and VBA?

texts:
Within the ESRI SDK resources there is an example of combining polygons by use of a geometry bag. The following code appears:

    'Define the spatial reference of the bag before adding geometries to it.
    geometryBag.SpatialReference = geoDataset.SpatialReference

That makes sense, based on some of the other documentation surrounding geometry bags, but when I tried it myself, I got a compile error saying that the method didn't exist. After looking at the geometry object model diagram, it appears as if the spatial reference property of a geometry bag is read-only... so what gives?


texts:
I'm developing an app with turn-by-turn navigation and I use OpenLayers consuming a tile server. However, I need to cache all the tiles of the route to be able to navigate offline. Is there a way to get all the tiles based on a LineString?

Currently just caching the tiles from the current position, but I need to get all the tiles in the route. I tried the following approach I found in the documentation, but as I understand it it is only caching only the seed current extent.

// start seeding the cache
    function startSeeding() {
        var layer = map.baseLayer,
            zoom = map.getZoom();
        seeding = {
            zoom: zoom,
            extent: map.getExtent(),
            center: map.getCenter(),
            cacheWriteActive: cacheWrite.active,
            buffer: layer.buffer,
            layer: layer
        };
        // make sure the next setCenter triggers a load
        map.zoomTo(zoom === layer.numZoomLevels-1 ? zoom - 1 : zoom + 1);
        // turn on cache writing
        cacheWrite.activate();
        // turn off cache reading
        cacheRead1.deactivate();
        cacheRead2.deactivate();

        layer.events.register("loadend", null, seed);

        // start seeding
        map.setCenter(seeding.center, zoom);
    }

function seed() {
        var layer = seeding.layer;
        var tileWidth = layer.tileSize.w;
        var nextZoom = map.getZoom() + 1;
        var extentWidth = seeding.extent.getWidth() / map.getResolutionForZoom(nextZoom);
        // adjust the layer's buffer size so we don't have to pan
        layer.buffer = Math.ceil((extentWidth / tileWidth - map.getSize().w / tileWidth) / 2);
        map.zoomIn();
        if (nextZoom === layer.numZoomLevels-1) {
            stopSeeding();
        }
    }

    // stop seeding (when done or when cache is full)
    function stopSeeding() {
        // we're done - restore previous settings
        seeding.layer.events.unregister("loadend", null, seed);
        seeding.layer.buffer = seeding.buffer;
        map.setCenter(seeding.center, seeding.zoom);
        if (!seeding.cacheWriteActive) {
            cacheWrite.deactivate();
        }
        if (read.checked) {
            setType();
        }
        seeding = false;
    }


texts: Dividing MultipPolygon using Polygon Divider

texts: How can I divide my shapefile polygon into smaller equal size area(20km2) using Polygon Divider, whenever I am trying to do a same size polygon shapefile generated?


texts: Python - Extract raster values at point locations

texts:
I'm a beginner with python. I would like to extract raster (.tif) values from a set of specific points that I have in a shapefile.
I would like to have the same result as the QGIS Point Sampling Tool plugin.
I can't install on Python libraries like qgis.core arcpy pygrass, but I have available GeoPandas GDAL Fiona Shapely rasterio and others.

How can I do this?


texts: QGIS QT .ui not displaying attribute table value

texts:
Following tutorial 6.3.6, Creating a New Form (the first introduction of Qt designer).

Seems very basic. Create a point layer with two attributes: name (text 80) and age (text 80). Then in QT designer create a .ui with name (QLineEdit) and age (QSpinBox)

I go to use the .ui in QGIS.

The form renders (but free space around the QLineEdit seems too large and very clunky).
Name - is populated from values from the table.
Age - is populated with 0, does not show the value from the table. (very odd)

Edits to name and age work fine, they are retained back in the table fine. Its just the initial display of age as 0 (not the age from the table) that's odd.

The autogenerated form works fine, both fields are populated with values from the table.

Seems odd that age defaults to 0 in QT? Perhaps this is by design with a spinbox default?

Brand new to QT designer (coming from a MSFT VBA environment for forms). I note there are 1500 pages of documentation on QT designer (O.K. that's pretty intimidating.) Wondering with ample experience in creating UI in MSFT how long does it take to become productive on the QT side? Are there any quick tips, posts, "got yas", tutorials, etc. to get productive on QT designer quickly?


texts: How to get co ordinates of points on mouse click in PyQGIS

texts:
I want to get x/y co-ordinates of point on mouse click for a point and need to use those co-ordinates to evaluate other value as per others module already define .

How can I get x y coordinates on mouse click in QGIS and print those in PyQGIS.


texts: ModelBuilder Iterator - how to change values and names respectively

texts:
I have a table (value.csv) with two fields (Expression, Name) that I would like to iterate through within a model.

I want the model to iterate through a series of SQL expressions (e.g. Value > 1, Value > 2, etc.). I think this part is working.

How do I get the iterator to rename each new shapefile according to the Name field in my table (i.e. Name1, Name2, etc.)

EDIT: I've added a screenshot for the set null parameters (which I am currently doing manually).

UPDATE: Added the Iterate Row Selection and Get Field Value to my model and after some tinkering it eventually worked. Here is the full model:


texts: How to create a custom Balloon for Google Earth?

texts:
I just came across one question related to google maps in stackoverflow.

Google Maps: How to create a custom InfoWindow?

I really like the popup that they mentioned in the question and the image is as follows.

Question: How to create or use same popup(InfoWindow) in google Earth.Is there any way to achieve this?

IMAGE:


texts: QGIS - Layers with CRS EPSG:3857 incompatible with Google street map

texts:
I have been having problems applying the different CRS systems to my layers. I have generate a variable distance buffer that I want to be circles around point on a Google Street Map of a city.

When I set the the CRS as 'EPSG:4326, WS 84' the circles appear as ovals:

When I set the CRS as 'EPSG:3857 Pseudo Mercator', while the shapes are the correct shape, they will not appear with the base map, and the base map is nowhere to be found.

This is confusing, given the the Google maps layer's CRS is 'EPSG:3857 Pseudo Mercator'!

I have read all the other questions/answers that I could find on this topic. I have tried saving it as a new layer (with different CRS) and changing the CRS of the layer before creating the buffer circles.

Edit: solution below


texts:
Is there a way to use a SELECT TOP (n) within ArcGIS desktop?  Is appears that since the ArcGIS query builder is effectively building the SQL statement after the WHERE clause, that this may not be possible.

Here is a typical SQL example of usage:

  SELECT TOP 3 * FROM Layer


texts: SaaS provider of overpass server for more that 10,000 request per day

texts: I am looking for a provider that runs overpass API that my company can use. We currently have about 10,000 request per day, but hope to grow in the future.


texts:
I have been searching for quality elevation data for the entire Scotland but I have not succeeded yet.

Do you know where to get the data in Raster or Vector formats?

Accuracy is very important as I intend to compare with other data.


texts: I am working through "Learn QGIS" for 3.4.  Under populating attribute data (pg 62) instruction say to use Field Calculator to populate area field of the attribute table.  Following instructions to the letter, I hit OK on the Field Calculator and, while the button appears to be active, nothing changes in the attribute field.  Is this a glitch?


texts:
I need to split the values in a field according to what they represent to be able to do a join on the data.

As seen above I have a ID field that needs to be split into the following:

The ID consist of 21 digits where the first 8 digits represents the Town Code

The next 8 digits represents the ERF number

The last 5 digits represents the Portion number.

In there a way to easily split these into their own fields as I need to join a excel sheet to Town number, ERF and Portion.

For example:

C04100030000194800000

Will be split into their own fields:

TOWN: C0410003

ERF: 00001948

PORTION: 00000

I assume that this will be done with arcpy as I have over 50 000 records?


texts: Why do I get different results running a script in PyWin and loading it to the Python shell in ArcGIS?

texts:
I've been trying to teach myself Python scripting for ArcGIS and have gotten frustrated on multiple occasions when I try to debug my script endlessly, only to later find out that it actually works properly when it is loaded into the Python shell in ArcGIS. In the past I have received a variety of errors when running the script in PythonWin, only to have the code work just fine when loaded into ArcGIS.

Why is there a discrepancy between the results of a script run in PythonWin and the Python shell in ArcGIS?

I am using:

Windows 8.1
Python 2.7.2
PythonWin (pywin 32 build 214)
ArcGIS 10.1 Advanced

The script is as follows:

#    Import arcpy, env
import arcpy
from arcpy import env

##   Enable file overwrite
env.overwriteOuput = True

#   Set workspace
env.workspace = "c:/documents/arcgis/arcgis tutorials/psu tutorials/lesson3practiceexercisea/washington.gdb"

#   Create new layer feature classes: Cities, ParkAndRide
try:
    arcpy.MakeFeatureLayer_management("CityBoundaries", "CityBoundariesLayer")
    arcpy.MakeFeatureLayer_management("ParkAndRide", "ParkAndRideLayer")
except:
    print "Could not make new layer."

I initially did not post the code because I've had this issue before with code of various types and I was thinking this was not a code error. This code works when loaded the the Python console in ArcGIS, but now when run from PythonWin.


texts: How to get all the lat/long coords that define a shapefile

texts:
I have a Shapefile with several shapes inside it.

Using MapWindow I have been able to determine that the shape I need is the one with shapeId 19.

As you can probably tell I have very little experience with GIS in general but I think I have come to the right place to ask for help.

What I need is to extract all the lat/long coordinates that determine the shape.

Is it possible to do this using MapWindow or do I need to use some other software?

If more information is needed to help, please comment and I will update as soon as possible.

Any help is appreciated as this is driving me crazy!


texts: Insert/Update/Delete table from another table changes in PostGIS

texts:
I'm trying to build a function that Insert/Update/Delete rows of a table based from another table changes. Unfortunately I could not continue, some problems occured:

DELETE isn't work correctly: Any linestring cannot deleted when trigger is enabled;
INSERT/UPDATE - More than 1 row is returned from  query - I insert the 'DISTINCT' to fix - I don't know why;
I would like to add point's manually (If possible).

The intention is create, update or delete a point for each line inserted, updated or deleted
when this line have a specific atribute called 'con';

These points are generated on the center of each line;

Point table - prec_peca_receive ;

Line table -  string

Here is the code:

CREATE OR REPLACE FUNCTION public.insert_peca_by_line()
        RETURNS trigger
        LANGUAGE plpgsql
        AS $func$

BEGIN
    CASE TG_OP
    WHEN 'INSERT' THEN
    INSERT INTO public.prec_point_receive (geom)

        select distinct st_lineinterpolatepoint(new.geom, 0.5)::geometry (point,31984)
        from public.string s
        where s.geom is not null and s.str_type = 'con';

    WHEN 'DELETE' THEN

        delete from public.prec_point_receive p
        where st_lineinterpolatepoint(old.geom, 0.5) = p.geom;

    WHEN 'UPDATE' THEN
    UPDATE public.prec_point_receive p

        set (geom) =
            (
            select distinct st_lineinterpolatepoint(new.geom, 0.5)::geometry (point,31984)
            from string s
            where p.geom is distinct from (st_lineinterpolatepoint(new.geom, 0.5))
            and st_lineinterpolatepoint(old.geom, 0.5) in
            (select p.geom
            from prec_point_receive)
            );
    ELSE
    --nothing

    END CASE;
    RETURN NEW;

END;
$func$;

CREATE TRIGGER tg_peca_by_line
    BEFORE INSERT OR DELETE OR UPDATE
    ON public.string
    FOR EACH ROW
    EXECUTE PROCEDURE public.insert_peca_by_line();

Could someone clarify?


texts: Clear or Disable Cache in ArcObjects Globe Control

texts: Is it possible to clear or disable the cache for individual layers on the AxGlobeControl? If so, how? I would like to be able to do this programatically through ArcObjects.


texts: Run command line shp2pgsql in C# function?

texts:
I'm developing an application web with ASP.NET/C# and I would in a button importing the shapefiles into PostGIS, and I'm trying with this code, this code I'm using it in cmd.exe in Windows it imported but in the application no but it can't import the shapefiles and it can't send me what'is the error, can you have idea?

if (strFileType == ".shp")
    {
    txtFilePath.SaveAs(Server.MapPath("~/UploadedShap/" + strFileName + strFileType));
    string strNewPath = Server.MapPath("~/Uploadedshap/" + strFileName + strFileType);

    // Connection String to shapefiles
    string[] strCmdText = {"cd ..","cd ..","cd program files","cd postgresql","cd 9.4","cd bin","shp2pgsql -s 3857 '"+ strNewPath+ "' public.COMGEO | psql -h localhost -d is -U postgres -p 5432 -w alami"};
    System.Diagnostics.Process.Start("CMD.exe",String.Join(" ",strCmdText));
    }


texts:
I'm an R user and have been using libraries like gistools and rgdal. I notice that ArcGIS appears to have features to add pie charts onto the maps themselves. (For example, if I have a map of the US, I can add a pie chart to each state). However, a quick search of gistools and rgdal documentation does not show any function to add such features.

Is there any easy way to add pie charts onto maps in R?


texts: pass string varialbe to placeholder for attribute in infoWindow (ArcGIS Javascript)

texts:
I am working on a Javascript-based webapp that uses an ArcGIS map service, which has several attribute columns containing URLs to photos.  The user will be able to open photos by clicking on a link in the infoWindow, accessed with an IdentifyTask.   I would like to use a variable in the placeholder.  My problem is I can't find a way to pass a string variable into the placeholder position.  Below is an example:

In my code, I have onclick=\"javascript:showPicture('${PrePic1}
where HTTP is the name of a column containing photo links.

I would like to use a variable, such as 'testString' so that I may vary the attribute in the placeholder (for instance, PrePic2, PrePic3 etc..)

I have tried setting a var to just 'PrePic1" and placing this in the brackets, and have tried '${PrePic1}' (with and without escape characters)..but these attempts didn't work..

any ideas? I hope I am explaining the problem well // Jason


texts: Can I use TIFF files for zonal statistics in QGIS or ARCGIS?

texts:
I am completing a remote sensing paper and am doing most of the work in PCI Geomatica, where I have 2 separate files. One from running a texture analysis algorithm containing 6 band layers from a satellite image. The other set from running an unsupervised classification on the same image with the original 6 band layers, as well as the texture analysis bands amalgamated to aid the classification, classification results which are two bands (input, output channels) giving me a total of 16 band layers. The classification itself is aggregated into 3 classes. My desire is to compare the correlation between the results of the classification (3 classes) against the results of the texture analysis quantitatively.

I am aware the GIS does not like pix files. I am considering converting both to TIFF files. My query is do I need to take any additional steps to run zonal statistics on these two image files correctly?


texts: What is elevation relative to in EPSG:2284?

texts:
I have some locations in EPSG:2284:

(3669486.63386,11293013.10592,709.08004),
(3669559.13811,11292972.72831,711.60055),
(3669639.51308,11292851.48264,712.22258),
(3669800.62304,11292949.38118,714.75766),

These are in US Survey Ft, however I don't understand what the elevation is relative to?

The reason this confused me if I have the following image showing these locations from the air and the EXIF data in the image reports an Absolute Altitude (which I assume is height above ellipsoid) of approximately 200m  - which would make it just under these points?


texts: Checking for number of shapefiles in directory using Python?

texts:
I want to test:

to see whether or not ONE AND ONLY ONE shapefile (which will be based on a filename wildcard) exists in a directory.
Add these shapefiles to a list
Use this list variable as an input to merge_managment function.

I already know how to walk through directories and pick out those filenames but not how to do the above.

I was thinking something like:

rootDir = "\homePath"
featureinput1 = "path"
newOnlyShpList = []
for dirPath, dirNames, fileNames in os.walk(rootDir):
    for file in fileNames:
       if fnmatch.fnmatch(file, "*" + "2012" + "*" + ".shp")== True ##AND IF ONE SHAPEFILE EXISTS BY ITSELF IN DIRECTORY:##
          ##append to newOnlyShpList
          arcpy.Merge_managment([featureinput1, newOnlyShpList], "OutputHere")
       elif
          ##ignore if greater than one shapefile

All I know about is the arcpy .exists function but all that does is return a boolean on whether any shapefile exists.  Maybe I need to check just for 1 .shp extension and it has nothing to do with arcpy?


texts: To filter gps co-ordinates that matches with the previous layer

texts:
I have a two .csv files with sets of gps co-ordinates in it. After plotting the gps co-ordinates from two files as separate layers, I need to filter a particular road which matches from the two layers and remove the remaining co-ordinates from the attribute table and save these co-ordinates back as a .csv file.

How to do this?


texts: Getting z coordinates from polygon list of Geodataframe?

texts: I am using Python to solve a problem with many polygons that could be really stressful with other software related to my field (mining), however i havent found a way to extract the Z coordinates from the Shapefile that I am working with.


texts: Displaying overlaping polygons and popups in the qgis2web plugin?

texts:
I'm total noobie when it comes to using qgis2web, CSS and HTML. The issue I'm having is displaying overlapping multipart polygons within the same shapefile (some of which are the same shape, but with different attribute data). In QGIS I have ordered the polygons to appear with the largest polygons at the bottom and the smallest towards the top. Some of the polygons do not appear at all. Where the polygons are the same shape it is difficult to tell if they have all come through, but only a single popup appears with one set of attribute data.

I have looked for answers on here, and I found this one: How can I get a QGIS2Web project to display info for multiple layers in a popup? but when I applied the code (which may have been in the wrong section), it stopped anything displaying at all from the qgis2web files.

I'm using QGIS 2.18.11 to create the shapefiles and qgis2web 2.22.0 to create the webmapping for wordpress.org.


texts: Exporting large map area in several smaller pages?

texts:
I have a large area at the scale 1:5000 which fits in a size larger than A0. I want to export it from ArcGIS into A3 pages at 1:5000. Below is in detail the exact area of what i am looking for and the position of the pages (which overlap) http://prntscr.com/cynt24 .

Is there a way to do this?


texts: We have a structured GIS from GeoServer version 2.15 (some layers loaded in shapefile format and others loaded from PostgreSQL). Therefore, we work with the GIS connected directly to the server and it works correctly, the problem arises when we connect from VPN, where the layers on the map are not shown (the layers loaded directly as shapefile and those that are in the base of PostgreSQL data). In this case, how can this situation be corrected and improved, from the point of view of the GeoServer configuration towards the VPN.


texts:
I am trying to install PDAL from github according to the following instructions:

Clone into a source directory: git clone https://github.com/PDAL/PDAL PDAL
Make a build directory: mkdir PDAL-build
Enter the build directory: cd PDAL-build
Run CMake to find dependencies: cmake ../PDAL
Or, see https://github.com/PDAL/PDAL/blob/master/cmake/examples/pramsey-config.sh
If dependencies are not found, manually set them: ccmake ../PDAL
Once CMake has found all dependencies, run the build: make all
And install the artifacts: make install

and I am kind of stuck. I don't really understand how the cmake prosses works. I have cmake and i do exactly what it says but I get a whole page of errors shown in the image. Anyone can help me or explain what I am doing wrong?


texts: Rasterize vector based on list values in Google Earth Engine

texts:
I would like to convert polygon (as featurecollection) that has the attribute 'Name' which contains string value of plantation types, namely 'Rubber', 'Mango', 'Coffee'. However, the .reduceToImage function won't accept string values. Thus, I had to do some workaround by deriving unique numeric id but now I'm stuck.

// Workflow: 1) Import asset, plantation_trial1
// print(plantations);
print(plantations.aggregate_array("Name"));

var plant_field = 'Name'
var all_text = 'what is'
// Function to build states list
var build_select = function(feature_collection, field_name) {
  var field_list = ee.Dictionary(feature_collection.aggregate_histogram(field_name))
                     .keys();
  return field_list.map(function(name) {
    return ee.Dictionary({'label': name})
  }).getInfo();
};
var plant_list = build_select(plantations, plant_field);

print(plant_list);

The 'plant_list' already contains the unique id for each 'Name', but no array or list function worked with the object. So basically, I am trying to have the script automate the extraction of the unique names from the featureCollection and assign each with a unique id. I am also aware that there may be some more efficient alternatives to achieve this, like distinct. Currently, here is what I got from print():

Can anyone give any idea on identifying what object that I have here and what to do with it, please?


texts:
I am trying to snap 2 polylines features as the picture below :

I would like to snap both end lines using tolerance distance ans only snap the nearest feature: from the red feature, choose the nearest line to the yellow one and snap it (extend it to the yellow end point).

I tried to use classic snapping and "new" snapping environment from ArcGIS editing session, but it works manually, and I would like to make this process for thousands of features...
That is why I also tried Snap Tool into Editing ToolBox, but this tool does not consider end points of the source feature, but feature vertices!

Does anyone has an answer which could save me from this problem please?


texts: How to use the IRelationalOperator interface to check if a point feature within a polygon feature?

texts:
I use the following code to check if a point feature is within a polygon feature

pFeaBlock is a polygon feature and Fea is a point feature

int blockID = -1;
IFeatureCursor pBlockCursor = this._blockLayer.Search(null, false);
IFeature pFeaBlock = pBlockCursor.NextFeature();
while (pFeaBlock != null)
{
    if (!pFea.Shape.IsEmpty)
    {
      IPolygon pPolygon = (IPolygon)pFea.Shape;
      IRelationalOperator pRelOperator = (IRelationalOperator)pPolygon;
      if (pRelOperator.Within(pFeaBlock.Shape) == true)
      {
          blockID = (int)pFeaBlock.get_Value(2);
      }
     }

     pFeaBlock = pBlockCursor.NextFeature();
 }

but I always get he COMException on

if (pRelOperator.Within(pFeaBlock.Shape) == true)

Is anyone familiar this interface?
Thanks in advance


texts: PostGIS select points where they intersect a polygon, but on some attribute values, apply buffer prior to intersect

texts:
I want to select some points based on a intersection with a polygon coming from an other database as a WKT geometry.
For the moment, everything is working fine:
SELECT p.name, p.geom, p.class
FROM point_table p
WHERE ST_Intersects(
        ST_Transform(p.geom, 2056),
        ST_Transform(ST_GeomFromText('POLYGON(((6.12 46.22,...,6.12 46.22)))', 4326), 2056)
      )
AND (...other conditions...);

But let's say the point have a 'class' attribute depicting the object they represent.
When class='city' I not only want to select the point which are inside the polygon, but also those which are in a 1000m buffer and for class='Town' I'd like to apply a 500m buffer.
I can do this in two different ways; either buffering the polygon itself, or the points, where the class attribute of the points matches city or town using the appropriate distance for each.
As the polygon is quite large, I prefer to buffer the points themselves.
I was thinking of something like:
SELECT p.name, p.geom, p.class
FROM point_table p
WHERE
  CASE
    WHEN p.class = 'City'
      THEN ST_Intersects(
             ST_Buffer(ST_Transform(p.geom, 2056), 1000),
             ST_Transform(ST_GeomFromText('POLYGON(((6.12 46.22,...,6.12 46.22)))', 4326), 2056)
           )
    WHEN p.class = 'Town'
      THEN ST_Intersects(
             ST_Buffer(ST_Transform(p.geom, 2056), 500),
             ST_Transform(ST_GeomFromText('POLYGON(((6.12 46.22,...,6.12 46.22)))', 4326), 2056)
           )
    ELSE ST_Intersects(
           ST_Transform(p.geom, 2056),
           ST_Transform(ST_GeomFromText('POLYGON(((6.12 46.22,...,6.12 46.22)))', 4326), 2056)
         )
  END
AND (...other conditions...);

I think this actually works, but I'm wondering if there is a better approach?
For example, one which would not repeat the massive polygon 3x.
I always want to "pre-compute" the geometry in my SELECT statement, using an alias, and then using this, but it's not interpreted in the rest of the query as if it doesn't exist at that moment, but only at the very end of the query:
SELECT p.name, p.geom, p.class,
CASE
  WHEN p.class = 'City'
    THEN ST_Buffer(ST_Transform(p.geom, 2056), 1000)
  WHEN p.class = 'Town'
    THEN ST_Buffer(ST_Transform(p.geom, 2056), 500)
  ELSE p.geom
END new_geom --<---- That would be a nice solution, but I cannot use 'new_geom' hereafter.
FROM point_table p
WHERE ST_Intersects(
        new_geom,
        ST_Transform(ST_GeomFromText('POLYGON(((6.12 46.22,...,6.12 46.22)))', 4326), 2056)
      )
AND (...other conditions...);


texts: How to define polygons overlapping each other?

texts:
I have a shapefile of landuse in different catchments. I want to calculate the percentage of each landuse category in each catchment using Tabulate Intersection. Before using this tool, I had dissolved the polygons in landuse shapefile using dissolve tool, dissolve field as the landuse category, checked create multipart feature, unchecked unsplit lines.

After using Tabulate intersection, I found out that the percentage of total landuse in some catchments are larger than 100% (around 103%). I think this means that polygons are overlapping each others in these catchment. As a result, I'd like to identify where polygons overlap each others to fix these errors.


texts: How can I control a custom parameter in QGIS WMS/WMTS/WFS webservices

texts:
GeoServer layers can be configured to be very flexible with parameters in at least these two ways:

viewparams pass through to the SQL view (link)
time/elevation/custom dimension (link and link)

A sample call would look something like: http://localhost:8080/geoserver/topp/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&STYLES&LAYERS=topp%3Atasmania_cities&SRS=EPSG%3A4326&WIDTH=768&HEIGHT=641&BBOX=145.9588623046875%2C-43.378143310546875%2C148.0682373046875%2C-41.617584228515625&time=2012-10-13T00:00:00.000Z&elevation=1233.47&viewparams=low:20000000 (note: time & elevation & viewparams=low)
In a JavaScript client like OpenLayers or Leaflet, it is easy to add these additional parameters to the WMS request.
In QGIS, it is possible to use the Time Manager (up to v3.12) or the Temporal Controller (after v3.12) to manipulate the time parameter (but the time parameter must be in the correct ISO8601 timestamp format, so can't be used for general parameters)

However, I cannot find another QGIS plugin that allows a user to manipulate any dimension (e.g. elevation or the low viewparam above) with a slider or similar control.
Can anyone point me in the right direction?
If not, is this something useful that I should consider implementing as a plugin?
I assume that I should be able to intercept the WMS/WMTS/WFS webservice requests when they are created, modify the parameters and forward them on.


texts: How to deactivate the Color Ramp in QGIS graduated vector styling?

texts:
I want to display points with different size and I DO NOT want to have a color ramp.
Do you know how I can do that?


texts: Select evenly distirbuted points of unevenly distributed set

texts:
I  have a task which looks pretty nontrivial yet pretty standard.

Let's say, I have a set of points, unevenly distributed on the map (i.e. in some regions there are a lot of points, while in the other regions points are sparse). The points can be sorted by some criterion K. I want to limit selection of point based on criterion K so that those points would be distributed relatively evenly. I.e. in the dense regions only really high on K points are shown, while in sparse regions whatever points there are are shown.

It seems to be a pretty standard task: e.g. showing cafes by their ranking or the most attractive tourists points of interest. So, I wonder, if there is a standard algorithm for this task (or maybe even a function in PostGIS) that I just not aware of?


texts:
I'm trying to create a barrier (rivers) that can't be crossed by the least cost path.

I've created a transition layer using "areas" which contains the transition values on land, and on rivers. However, the barrier has a 'buffer' around it, with increasing values, instead of being a discrete barrier.

Cost raster and transition layer

As can be seen, there's a 'buffer' around the barrier in the transition layer.

Any ideas on how to create a discrete boundary? preferably with the barrier as a value 1000, with everywhere else a value of 1.

code:
library(gdistance)
    ## Create cost surface where "land" exists in the middle
    cost <- raster(nrow=100, ncol=100,
               xmn=0, xmx=100, ymn=0, ymx=100, crs="+proj=utm")
    cost[] <- 10
    cost[cellFromRowColCombine(cost, 50:55,20:80)] <- 1000
    costf <- asFactor(cost < 100)
    ## Produce transition matrices, and correct because 8 directions
    trCost <- transition(costf, "areas", directions=16)
    trCost1 <- geoCorrection(trCost[[2]], type="c")
    plot(cost)
    plot(raster(trCost[[1]]))


texts: Hydrologic analysis with sextante-saga in qgis

texts:
(rough translation): Hi, I want to do hydrologic analysis of a DEM with sextante-saga in qgis. I found the plugin and the needed algorithms in SAGA, enabled them in the sextante toolbox, but when I run them, they error out saying SAGA is not installed (not where sextante is looking for them).

My question is therefore where to find these SAGA binaries? Thanks

  Salve, dal nome penso che sia italiano come me. Spero possa darmi una
  mano... Voglio fare l'analisi idrologica da dem su qgis (1.8) con
  sextante. Ho scaricato il plugin e messo in
  C:\OSGeo4W\apps\qgis\python\plugins .. l'algoritmo di flow direction
  (catchment area) e gli altri algoritmi sono in SAGA, che ho attivato
  dalla toolbox di sextante tramite interfaccia qgisdesktop, quando
  lancio l'alforitmo mi viene detto che devo prima installare saga, ho
  capito che devo scaricare a parte SAGA binaries e metterlo sotto la
  folder di sextante... Può indicarmi dove posso scaricare
  quest'applicazione? In un suo post passato ho visto che lo ha chiamato
  msw SAGA package! Grazie INFINITE per la dritta! Cordialmente laura guarnieri


texts: Using Copy Raster in ArcPy gives different result to tool dialog?

texts:
I am trying to convert 16 unsigned bit tifs to 8 unsigned bit BMPs. When I used the Export Raster in ArcGIS, it works well. But when I used arcpy

arcpy.CopyRaster_management(each,
                            "E:/HutDetection/17.bmp","DEFAULTS",
                            "","","","","8_BIT_UNSIGNED","","","BMP","")

the output images is too bright and its pixel value is not the same as the one got from ArcGIS.
I have changed the parameters a lot of times but I cant get the result I want.

What is the problem?

result from arcpy:

result from arcgis: