texts: Divide multilineZM into parts of same length

texts:
I have a postgis table of multilineZM features. I need to divide them into lines of same length (eg 1000 meters). I am new to QGIS and PostGIS (I usually use ArcGIS and ArcObjects). I found in the documentation of QGIS the function QgsGeometry::locateBetweenMeasures but it seems doesn't work with M coord. There is a similar function in PostGIS: ST_LocateBetween that works with M but I don't know how to use it with python QGIS. I tried to define a View in postgres, but DBManager says it's not a valid layer (the same query in postgres sql editor works).

The query is:

SELECT
    myTable.gid,
    st_locatebetween(myTable.geom, 0::double precision, 500::double precision)
AS
    myNewTable
FROM
    myTable
WHERE
    myTable.gid = 19633;


texts: Can TileMill be configured to produce vector tiles, given a mapnik that supports the format?

texts:
Out of the box, TileMill 0.10.2 can import shapefiles from OpenStreetMaps (OSM), can style layers with CartoCSS, and can output an mbtiles file containing raster tiles (png). TileMill uses mapnik as part of its processing pipeline to go from map-data to rasterized map images (the output).

As of now, April 2017, it is possible to modify Mapnik with a library/plugin to produce Mapnik Vector Tiles (MVT) files instead of raster images. The project is here: mapnik-vector-tile (github).

Can we substitute or augment the "Vanilla" mapnik used in TileMill with the newer version?

Note: A different-but-related question was asked two years ago. It asks if MapNik can output vector tiles. Yes it can, if that library is installed. What I'm asking here, is "given that one can force mapnik to spit out vector tiles, how can I make TileMill call into the modified library?". Is there a way to "shim" mapnik to output MVT tiles?


texts: Problems with weighted overlay on QGIS for Mac

texts:
I am a beginner on GIS seeking to finish a suitability analysis. For this I was planning to use 'vector to raster' tools as well as a 'weighted overlay analysis' to find the most suitable areas.

I was able to rasterize a vector layer in QGIS 2.14 Essen, using raster > conversion > rasterize.

However I can't find a way to do the weighted overlay analysis. All threads I've read says that SAGA plugin should have such a tool, but I cannot find SAGA in neither QGIS versions 2.14 or 2.18.

I have Mac OSx version Sierra. I have searched for weighted overlay tools on both QGIS 2.14 Essen and 2.18 Las Palmas without luck. I don't have time to downgrade my MacOS and hope the problem is not with Sierra compatibility.

How can I get the SAGA plugin to work? Or in case its impossible, what other solutions are there?


texts: Enabled Level FINE logger GeoTools/GeoServer

texts:
I would like to enable my logger in level FINE with GeoTools.
I changed my DEFAULT_LOGGING.properties to :
log4j.category.org.geotools=DEBUG
log4j.category.org.geotools.factory=DEBUG
log4j.category.org.geoserver=DEBUG
log4j.category.org.vfny.geoserver=DEBUG

But still when I debug geotools lib, LOGGER.isLoggable(Level.FINE) always return false and I don't know, why ?


texts: Postgis query on 3D cube does not use index

texts:
I am testing Postgis on the follow query:

SELECT * FROM particles
WHERE position  &&&
'LINESTRING(0.08537 0.0738 0.0145,0.0953 0.0837 0.0345)';

that returns about 100k rows from a data set of about 100M. The schema of the table is

Table "public.particles"
  Column  |       Type       | Modifiers
----------+------------------+-----------
 partid   | integer          | not null
 time     | double precision | not null
 position | geometry         |
Indexes:
    "particles_pkey" PRIMARY KEY, btree (partid, "time")
    "particles_the_geom_gist" gist ("position")
Check constraints:
    "enforce_dims_the_geom" CHECK (st_ndims("position") = 3)

Even tought I have an index on the attribute position, the query is always executed  with a full table scan and  so it is  particularly slow.

I also tried

set enable_seqscan = off;

but it did not change anything.

Using different query like,

SELECT count(*) FROM particles WHERE ST_3DDWithin(position, 'SRID=4326;POINT(0.08537 0.0738 0.0145)',0.01);

which retrieves objects within a give distance from a point, Postgres uses the index, why it doesn't for the other  query?


texts: How to conduct a query in another scene layer in a webscene?

texts:
I have the following situation,

I have imported a web scene into a script. The webscene consist of two scene layers and during most of the script I am only working with one scene layer. However almost at the end of my script, whenever I want to click on an anchor element in a list in the sidebar which I have created with HTML, I would like to zoom in on an object in another scene layer.

The anchor element does contain a dataset.gebouwid and I include this information in a new Query. However now I am completely lost in which way I can use this query to query through another sceneLayer.

I have looked at the sandbox of Query SceneLayerView, however in this specific case throughout the whole script there appears only one scene Layer. So I was wondering what in case when you want to query through another scene Layer.

 for (i = 0; i < uniqueItems.length; i++) {
                                var value2 = uniqueItems[i];
                                var subtext = value2.substring(value2.length - 9, value2.length);
                                //   console.log(subtext);

                                var entry = document.createElement('li');
                                var link = document.createElement('a');
                                link.textContent = value2;
                                // link.href=value2;
                                link.href = '#';
                                link.dataset.gebouwid = subtext;
                                // console.log(subtext);
                                link.onclick = function (event) {
                                    var clickedElement = event.target;
                                    var element = clickedElement.dataset.gebouwid;

                                    var query = new Query({
                                        where: "lokaalid = '" + element + "'",
                                        returnGeometry: true,
                                        outFields: ["*"]
                                    });


texts: Qgis project in Lizmap not showing layers, only base layers

texts:
I am following the lizmap tutorial for publish a QGIS project in lizmap but I only see the base layers. Checking and unchecking the layers don't show them. The demo Montpellier project is working well but not my project. Using lizmap 3.1.2 , QGIS 2.18, QGIS Server 2.18.10 on Ubuntu 16.04

This is my QGIS project

Any suggestions?


texts:
Here is the case:
Im creating a new marker on map:
centeredMarker = L.marker(map.layerPointToLatLng([mapCenterXPx, mapCenterYPx]), markerOptions).addTo(map);

And then it should center and stay fixed on view center of map when map dragging by user. Here is my implementation:
map.addEventListener('move', toCenterMarker)

where toCenterMarker function is
function toCenterMarker() {
  let mapCenterPx = map.latLngToLayerPoint([map.getCenter().lat, map.getCenter().lng]);
  let mapCenterXPx = mapCenterPx.x;
  let mapCenterYPx = mapCenterPx.y - 80;
  let mapCenter = map.layerPointToLatLng([mapCenterXPx, mapCenterYPx]);
  centeredMarker.setLatLng(mapCenter);
}

So the problem is when user dragging map, marker fixed on center and thats feels fine, but user can see that this marker have a little shakes on moving. Is there any solution of this problem?


texts: Why changes made in QtDeisgner form not shown in QGIS plugin gui?

texts:
I have created a plugin in QGIS 2.10.1. After making some modifications in the ui file in Qt Designer and reloaded the plugin or restarted QGIS or restarted the pc these are not incorporated. I cannot understand why, since till now this has been working. The only differences with other times is that before the modifications I had made a copy of the the plugin folder (copied folder), pasted it on another folder on my computer, then erased the initial plugin folder from .qgis2\python\plugins\ folder and re-pasted the copied folder inside the .qgis2\python\plugins\ folder. Due to my limited relevant knoweledge, I thought to compile again the resources file, but again no result. Also, I applied a make deploy (when trying to change the icon) so maybe this has to do with this? That' s all the actions I have taken. So, any idea why changes made in the QtDeisgner form are not shown in QGIS plugin gui?

But if the folder that conatins the modified ui, is copied and pasted on another folder in my computer, erased from the .qgis2\python\plugins\ folder in wich the folder that conatins the modified ui is pasted, when running again the QGIS I can see the changes. But again if I make new changes nor plugin reload neither restarting QGIS is working. I have to follow again the procedure I just described.

thanx


texts: Downloadable pdf format of GRASSGIS 7.0 manual

texts:
Can anyone help with a link or something that could enable the download of a pdf version of the grassgis manual? I am having a hard time trying to convert the html version into printable pdf.
Nkeze


texts: Need Image with Sweref99_1800 projection/coordinate system

texts:
I need third party image tiles (probably world imagery) under sweref99_1800 projection and coordinate system (EPSG 3008). I found this link http://www.naturalearthdata.com/downloads/10m-raster-data/10m-cross-blend-hypso/
for downloading world imagery but it is in spherical mercator projection (EPSG 4326). Is there any link from where we can download world imagery in sweref99_1800 format.

OR

is it possible to change the projection and coordinate system of image tiles (world imagery) from spherical mercator(EPSG 4326) to swereff99_1800 format (EPSG 3008) ?


texts: What prevents GRASS and SAGA integration into QGIS when users aren't administrators?

texts:
Windows version of QGIS 2.0 using the stand alone installer. I have been able to tie up the functionality using the Options and Configurations settings in the Processing menu.

This works fine in a standard PC and logged in as an "administrator". However something is preventing it from working with the same QGIS setup conditions when not logged in as administrator.

Can anyone shed any light on what might be happening here - especially what process the Processing tools go through when a non-QGIS module is activated. Non-administrators cannot write to the QGIS folder, nor can they put executable files into temporary folders. Could it be one of these issues?


texts:
I am trying to write any error messages that arise in my script and record it in a text file. When I check my errorLog, the error messages are not recorded. Below is my code: Note the parameter 'filePath' is the file path to the file that produced the error.

except:
    addError = arcpy.GetMessages(2)
    with open(errorLog, 'a') as errorMsg:
        errorMsg.write('%s, %s\n' % (filePath,addError))

Can anyone offer insight?


texts: Is there a Small Windows install of QGIS?

texts: I am looking to find a faster alternative to ArcReader.  I think QGIS might be faster but I was wondering if there was a smaller install without all the GRASS etc for windows.  I just want to get a fast map out to my users in the field. The ability to Zoom, Pan, Identify and perhaps search for an address would be nice.


texts: How to automatically adding z dimension in vertex editor from field attributes

texts: I have geo-coordinates x y and z value, when I created an 3D shapefile, how can I automatically make the z value into Vertex Editors' z dimension??


texts: Selecting the MAX value from a feature using ArcObjects?

texts:
I have a feature (point) with a field called 'Sequence'. The values in this column range from 1 to n. I can use IQueryFilter to get the first value, 1.

                    Dim pQry As IQueryFilter = New QueryFilter
                    pQry.WhereClause = """Seq"" = 1"
                    pCursor = pFC.Search(pQry, True)
                    If TypeOf pCursor Is IFeatureCursor Then
                        pFeat = pCursor.NextFeature
                        If TypeOf pFeat Is IFeature Then
                            fStartElevation = pFeat.Value(nElevationFld)
                        End If
                    End If

I cannot figure out how to use this method to select the MAX value, since it varies across features. Is it possible? I'm thinking I just need the correct syntax. I have tried a few variations of pQry.WhereClause = """MAX""Seq""" after doing some forum reading regarding aggregate functions in ArcGIS, but I haven't found a solution that works.

If this isn't possible (I figured it should operate similarly to SQL), is there another preferable method for performing this task? I am aware of Select By Attribute, but the syntax for this is tricky (examples online do not work in Arc10) and still requires the attribute table to be read.

EDIT: Here is a working solution using @George's suggestion based on IDataStatistics. I may just turn this into a small function that takes in a Featureclass and field name for the sake of re-usability.

                    pCursor = pFC.Search(Nothing, False)
                    If TypeOf pCursor Is IFeatureCursor Then
                        Dim pData As IDataStatistics = New DataStatisticsClass
                        pData.Field = "Sequence"
                        pData.Cursor = pCursor

                        Dim pStatResults As IStatisticsResults = pData.Statistics
                        Dim intMaxSequence As Integer = pStatResults.Maximum
                    End If


texts: Converting string column to datetime column in QGIS

texts:
I have an attribute column that has a string datetime and I want to convert that to a date column. I first attempted to simply add the time column in an expression, and I only got the year/month/date, but I wanted the time as well. So I first converted it the string to a new column with format_date(), to clean up the format, and then I went ahead and used to_datetime() to convert the string into the date column. It shows up fine in the Expression dialogue but then doesn't add the hour/minute/seconds.
What am I overlooking?
Using
QGIS version
3.16.1-Hannover
on Windows


texts: Calculating minimum distance between points and polygons in QGIS?

texts:
I have two vector layers: one point and one polygon and I want to calculate the minimum distance of each point from the polygons (ie the distance from the closest point of the closest polygon).

Is there any available plugin in QGIS for doing this?


texts: set wrong GPS unit in "degree" and want to change to "meter" in GIS while editing

texts:
I want to edit the GPS result in ArcGIS, but I found that there is a problem.
I set the wrong unit while I used the GPS outside , and the unit is "degree".
Although the coordinate system is correct, I still cannot change the unit to meter.

Is there any method that I can edit the unit to meter in GIS ?  (or QGIS?)
Thank you very much!


texts: Request WFS from GeoServer with Jquery AJAX using INTERSECT cql_filter error

texts:
I am having difficulty troubleshooting this issue and would appreciate some ideas. Thanks in advance!

I am trying to make an AJAX request from a Leaflet map with Jquery to GeoServer to return a WFS selection using the INTERSECT sql_filter param. I have it working fine using the BBOX sql_filter.

Successful request using BBOX filter

   var geojsonLayer = new L.GeoJSON();
   var currSelectionLayer = 'DRCOGPUB:rea_drcog_jurisdiction_county';
   var sw_lat = 39.93290692296977;
   var sw_long = -105.39459228515624;
   var ne_lat = 40.094882122321174;
   var ne_long = -105.150146484375;

   function loadGeoJson(data) {
            console.log(data);
            geojsonLayer.addData(data);
        }
   var geoJsonUrl = "http://gis.drcog.org/geoserver/DRCOGPUB/ows?service=WFS&version=1.0.0&request=GetFeature&typeName="+ currSelectionLayer +"&srsName=EPSG:4326&maxFeatures=10000&outputFormat=json&format_options=callback:getJson&cql_filter=BBOX(the_geom," + ne_long + "," + ne_lat + "," + sw_long + "," + sw_lat + ")";

   $.ajax({
            url: geoJsonUrl,
            dataType: 'jsonp',
            jsonpCallback: 'getJson',
            contentType: 'application/json',
            success: loadGeoJson
        });
   map.addLayer(geojsonLayer);

However, when I apply an INTERSECT filter I get this error:

  Resource interpreted as Script but transferred with MIME type text/xml

Here is my code for trying to get the INTERSECT filter to work:

Un-successful request using INTERSECT filter

        var geojsonLayer = new L.GeoJSON();
        var currSelectionLayer = 'DRCOGPUB:rea_drcog_jurisdiction_county';

        function loadGeoJson(data) {
            console.log(data);
            geojsonLayer.addData(data);
        }

        var geoJsonUrl = "http://gis.drcog.org/geoserver/DRCOGPUB/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=" + currSelectionLayer + "&srsName=EPSG:4326&maxFeatures=10000&outputFormat=json&format_options=callback:loadGeoJson&cql_filter=INTERSECT(the_geom,%20POINT%20(-105.36369323730467%2040.10013461308659)";
        $.ajax({
            url: geoJsonUrl,
            dataType: 'jsonp',
            jsonpCallback: 'getJson',
            contentType: 'application/json',
            success: loadGeoJson
        });
        map.addLayer(geojsonLayer);

I have tried many variations on the AJAX request and tweaking my Apache config but have had no luck.  Any ideas would be appreciated. Please let me know if you need more info.

UPDATE:

I enabled the mime_module (mod_mime.so) in Apache and added the content type with this config code:

AddType text/javascript .js

UPDATE - Final Fixed Code:

        var geojsonLayer = new L.GeoJSON();
        var currSelectionLayer = 'DRCOGPUB:rea_drcog_jurisdiction_county';

        function loadGeoJson(data) {
            console.log(data);
            geojsonLayer.addData(data);
        }

        // using INTERSECT strategy
        var geoJsonUrl = "http://gis.drcog.org/geoserver/DRCOGPUB/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=" + currSelectionLayer + "&srsName=EPSG:4326&maxFeatures=10000&outputFormat=json&format_options=callback:getJson&cql_filter=INTERSECTS(the_geom,%20POINT(-104.986145%2039.733066))";
        $.ajax({
            url: geoJsonUrl,
            dataType: 'jsonp',
            jsonpCallback: 'getJson',
            contentType: 'application/json',
            success: loadGeoJson
        });
        map.addLayer(geojsonLayer);


texts: Joining attributes by location/intersect works for all but one polygon in QGIS?

texts:
I used the Join attributes by location and Intersect tools in QGIS 3.4 for my points layer to add a column with name of the intersected area into the attribute table of points. It works perfectly for 5 out of my 6 areas - they are all in one shp with the same characteristics. I tried to export it separately, save it differently, but this one area never works.

Any ideas what could I try?

See printscreen - I selected the result I get after I use the tool and the area for which it doesn't work - the points intersected are not in the result layer.


texts: Download images from GEE directly to the computer

texts:
I have few codes that retreive LANDSAT images, calculate NDVI and then create rasters and then export them to my drive. Mnay times if the raster is too heavy it breaks it into few images.
My question is if there is any way to export images directly to my computer?

This is an example of code which I would like to downliad its' products to the compuer but not to my drive:

var geometry=plots;
/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000)
  .copyProperties(image, ['system:time_start']);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2019-01-01', '2019-09-12')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var clippedCol=dataset.map(function(im){
   return im.clip(geometry);
});

//test if clipping the image collection worked
Map.centerObject(geometry,14);
Map.addLayer(clippedCol.median(), rgbVis, 'RGB');

// Get the number of images.
var count = dataset.size();
print('Count: ',count);
// print(clippedCol);//here I get the error messege "collection query aborted after accumulation over 5000 elements
// print(dataset,'dataset');//the same error here

//function to calculate NDVI
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(ndvi);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

// // Test the addNDVI function on a single image.
// var ndvi1 = withNDVI.select('NDVI').mean();

var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

//get NDVI images

var listOfImages =(withNDVI.toList(withNDVI.size()));

var listOfNumbers = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,23,24,25,26,27,28,29,30,31,32,32,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,];

for (var i in listOfNumbers) {
  var image = ee.Image(listOfImages.get(listOfNumbers[i]));
  var toexport=image.visualize(colorizedVis).addBands(image);

  // export
  Map.addLayer(image, colorizedVis, i);
  Export.image.toDrive({
  image: toexport.toFloat(),
  description: i,
  scale:2.5,
  crs:'EPSG:4326',
  maxPixels:1310361348,
  region:geometry.geometry().bounds()

});

 }

//Long-Term Time Series

print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.mean(), 30));


texts: Manually editing raster to improve classification using ArcGIS Desktop?

texts:
I have a land cover classification derived from Landsat TM imagery - I created this classification using ENVI. I would like to edit individual and groups of cells in some areas of the classification in order to improve the accuracy of the output by comparing it with aerial imagery.

Is there a way, using ArcGIS Desktop, to manually edit large groups of cells by selecting them and assigning them to a new or existing class?

I am aware of merging methods based on class statistics, but the area I'm classifying has a high amount of spectral confusion due to thin cloud cover - it seems like a manual edit is the only way to go.


texts: Fortifying re-projected shapefile for ggplot2 in R

texts:
I have a shapefile of US county boundaries which I got from the NHGIS. I can load the file using the rgdal package and see information about it:

library(rgdal)
shp_1860 <- readOGR("shp", "US_county_1860", stringsAsFactors = F)
summary(shp_1860)

# Output:
Object of class SpatialPolygonsDataFrame
Coordinates:
       min     max
x -2356112 2258225
y -1337373 1565782
Is projected: TRUE
proj4string :
[+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0
+datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0]

You'll notice that the projection is a little unusual. But what I really don't understand is the latitude and longitude: what do those values mean? Why aren't they latitudes and longitudes in the vicinity of the United States?

I can fortify this shapefile for use in ggplot2, and plot it without a problem.

library(ggplot2)
shp_1860_df <- fortify(shp_1860, region="GISJOIN")
p <- ggplot() +
  geom_map(map = shp_1860_df, data= shp_1860_df,
           color = "black", fill = "white", size = 0.1,
           aes(map_id = id)) +
  expand_limits(x = shp_1860_df$long, y = shp_1860_df$lat) +
  coord_equal()
print(p)

I get a map just fine with that code. If I replaced coord_equal() with the more logical coord_map(), R and my computer always crashes because memory use skyrockets. I presume this is because of the weird coordinates.

I thought I would re-project the shapefile in order to plot it.

shp_1860_proj <- spTransform(shp_1860, CRS("+proj=longlat +datum=WGS84"))
summary(shp_1860_proj)

# Output
Object of class SpatialPolygonsDataFrame
Coordinates:
        min       max
x -124.7625 -66.94993
y   24.5211  49.38436
Is projected: FALSE
proj4string :
[+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]

Notice that the latitude and longitude now have values that make sense. I can plot this new shapefile using base R without a problem, e.g., plot(shp_1860_proj). But when I try to fortify it, I get an error:

shp_1860_proj_df <- fortify(shp_1860_proj, region="GISJOIN")
Loading required package: rgeos
rgeos version: 0.3-1, (SVN revision 413M)
 GEOS runtime version: 3.3.3-CAPI-1.7.4
 Polygon checking: TRUE

Error in createPolygonsComment(p) :
  rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 6

What is the right way to go about re-projecting and then fortifying this shapefile? Alternatively, can I re-project it after I've fortified it?


texts: Points in Polygon results in PNTCNT=NULL

texts:
I am trying to use the QGIS (2.6.1) vector analysis tool "Points in Polygon". The points originally came from a KML file and the polygons I created new. I initially had a warning about different CRS but resolved that by saving the files as shape files with ESPG:4326 WGS84 CRS. Now the tool runs fine but the PNTCNT field in the resulting shapefile is all NULLs.

Am I missing something? What should I check for?


texts: Using Pyqgis 3 to select features from attribute table by expressions

texts:
I have a shp with the following attribute table:

The column Line_ID represents each line (line 0, line 1, line 2 ... ).

The column ID are the points id (point 1, point 2 … ).

The column Declividad is the slope.

I need to find a specific slope break for each line.
So, the way I’m thinking to do this is to check if the previous point has slope value less than 8 and if the posterior point has slope value more than 8. So the script should select the posterior point.

For example: line 0 / point 1 = slope 2 ;line 0 /  point 2= slope 1;line 0 /  point 3= slope 4; line 0 / point 4= slope 10.
So the script, in this case, should select point 4, because it has slope value more than 8 and its previous point (point 3) has slope value less than 8.

How could I make this using Pyqgis?


texts: Unix date displayed for ESRI Collector date input?

texts: Trying to set up a date input for a Collector project on iOS, but the app wants dates in 13 digits (I believe Unix or Java time format?). How can I get the input to accept something human-friendly (i.e. m/d/y would be phenomenal).


texts: Removing noise selectively from DEM using QGIS?

texts:
I am looking for a processing tool that can remove noise selectively from a DEM (Mostly SRTM). My current workflow includes "Fill sinks" and "r.neighbours". The problem is that the tool is blind to places of abrupt slope changes and smooths the surface uniformly. I would like to apply a different amount of smoothing to areas of different slope values.

How can I do that?


texts: Can two polygon layers be merged in QGIS on a common id?

texts:
Can two polygon layers be merged in QGIS on the basis of a common id?

Update:

Ok , the dissolve tool works fine. But for instance lets say I have 10 polygons. 8 of them are together , so when dissolved they become one big polygon. The other two polygons are somewhere away and I dont want them to dissolve with these other 8. So how do you go about it.


texts: That is possible on a web map designed by R & leaflet give the opportunity to the user to draw lines, polygons, choose the colors (like a map editor)? if yes what's the package that allows me to do that?


texts: Join OS Building Heights and MasterMap Topography

texts:
I am trying to join the OS Building Heights CSV file with the OS MasterMap Topography Layer. They have both been paid / they are NOT the free OS data and so definitely include the data I need. This is in order to create a building height plan of my study area.

Having googled it numerous times, I am under the impression that I need to join the TOID field in the OS MasterMap Topo layer with the equivalent number in, what is called, 'field 1' of the building heights CSV file.

Problem is... when they join, all the 'new' fields which have been joined all say 'NULL'. An example of the number in 'field 1' is osgb1000002107214922 and an example of a number in TOID is 0001000034928367. So the Building Heights CSV file displays 'osgb' before the number and the OS Topo does not.

There is far too much data to go through and add in/remove the 'osgb', so does anyone have any other solutions or know where on earth I am going wrong?

Or could explain to me how to join these datasets in QGIS?


texts: Reprojecting sources with MapProxy yields distortion

texts:
I have configured MapProxy to aggregate imagery from a number of Internet sources so that I can dynamically overlay them by adjusting the layers=<x> query that I pass to MapProxy's WMS service. This all seems to be working properly.

Here's the problem: some of the sources are TMS servers, which use the spherical Mercator (EPSG:900913) projection. Some of the others are WMS servers that only serve data in the EPSG:4326 (WGS84) projection. MapProxy is nice enough to handle this for me, caching the images in their native projection and then reprojecting them to the target projection when I request them. However, there is a great deal of distortion in the resulting images near the -180/+180 degree longitude line. Here's an example:

The source image is from MapQuest's Open Aerial service, which is in global Mercator. The above is MapProxy's reprojection of the source to EPSG:4326. For most of the globe, it looks great, but the rightmost part of the image is obviously way off. If I request it in its native projection, then I get a nice image:

I see very similar results if I take an EPSG:4326 source (Landsat images from this WMS source) and request it from MapProxy using EPSG:900913:

Am I missing something fundamental here? Is this to be expected?


texts: ArcGIS for Server 10.2.1 cache scale greyed out in Service Editor and not all scales available in Manage Map Server Cache Tiles

texts:
I am using ArcInfo desktop 10.2.1, with an admin connection to ArcGIS for Server Standard 10.2.1

I have used the "Create Map Server Cache" tool in ArcToolbox to prepare all the underlying structure for the cache of an existing un-cached basemap service.

I've used a predefined tiling scheme, which tells it to generate cache tiles at 3 scales - 1:8000, 1:5000, 1:2000.

I've set the tile cache format the "PNG" and the Storage format to "Compact".

This tool runs without any problems.

The problem comes when I come to actually generate the cache using the "Manage Map Server Cache Tiles" tool from the toolbox.

Only 2 of my 3 scales are visible in the "Scales" box.

In Catalog, if I right click on my service > Service Properties > Caching, the slider shows the full extent of the 3 scales I want to cache at.

If I click on the Advanced Settings tab, I can see the 3 scales, however one of them is greyed out.  This is the one that's missing from the "Managae Map Server Cache Tiles" tool.

Is this a bug in the software? Or is there something I have missed?
I've also tried it with ArcInfo 10.1 and I get the same problem.

Is there a limit to the number of tiles a single cache scale can hold? Or a size limit?

If this is a bug, is there a workaround?


texts: Change Symbology with different Column Value

texts:
I am using ArcMap 10.0.

I am trying to change Cadastre colours when the values of the box EXTRA is set to 0,1 or 2.
So basically, I need the default values set to 0 (already completed).
Properties that have to be invested have a value set to 1 (completed).
Properties with further information needed value set to 2 (completed onsite).

To make it user friendly I would like to set up the symbology colour to change the polygons to different colours depending on the EXTRA value.

0 = clear
1 = green
2 = red

I have tried both Quantities and Multiple Attributes but neither get the result I am after.

Using the 'Graduated Colors', changing Fields Value to EXTRA only gives one symbol when I would want three? If I have three then I will be able to change the range and respective colours.


texts: Preparing CSV files for use in ArcGIS Desktop?

texts:
How do I prepare CSV files for use in ArcGIS Desktop.

I ask because I have some troubles using CSV files because ArcGIS attributes wrong field types to my columns and also misinterprets special characters such as á or ê.

I have read in the Esri forum that there is a so-called schema.ini file that defines somehow the field types e.g "Col22=V002 Text" see here http://forums.esri.com/Thread.asp?c=93&f=1149&t=64464

That's kind of funny because I have often seen these .ini files on my disc but never actually wondered what they are good for. It is kind of weird that Excel stores such metadata in an extra file since other programs like R don't do so.

I already tried to manipulate this .ini file with little success since I didn't find out how to apply for example "string" type. There are some information on MS sites, see here: http://msdn.microsoft.com/en-us/library/windows/desktop/ms709353%28v=vs.85%29.aspx but I couldn't find a solution.

Also I didn't really like the idea to work with this .ini file because it is quite a bunch of work to define and type all the fieldnames when I have let's say 50 columns. And these .ini files might get lost, etc.


texts:
Im running my script in ArcGIS 10.2 and the dialog window says completed, but the empty fields in the 'SoilsType' row of the 'mSoils' Feature Class aren't getting populated with the values from the 'SoilsType' field from the 'mLanduse' Feature Class.

cur = arcpy.da.SearchCursor(mLanduse, lField)
for row in cur:
 cur2 = arcpy.da.UpdateCursor(mSoils, sField)
 for row2 in cur2:
     if row2[0] == '' or row2[0] is None:
         row2[0] = lField
         cur2.updateRow(row2)


texts: How to get the tripoint of three or more polygons in PostGIS?

texts: Assume that I have a column containing polygons. What is the SQL statement to get the list of points that three or more polygons touch with each other?


texts: How can I Export transparency and text well when using qgis2web in QGIS 3.8?

texts:
Hi, I am an amateur who is making a historical map.
As you can see in the picture, I have a problem using qgis2web.

The first picture is a normal scene output from QGIS version 3.8.

The second picture was taken from qgis2web with a Mapbox. Transparency is nice, but the text in the labels is terrible.

The third picture was taken from qgis2web with openlayers. Transparency, text in labels, it's all a mess.

I tested this site a year ago, but it was fine at the time. What should I do ? It's too hard...


texts: How to create image pyramid for GeoServer programmatically

texts:
I try to create a tiled layer for my OpenLayers3 application. For the creation part I plan to use GDAL gdal_retile.py script. I've already tried it and it seems to work as I expected it to. Now I want to provide the resulting tiles as a layer to my OpenLayers application. Indeed (I guess like all programmers), I do not want to use GeoServer UI for this part and want to do it compeletly programmatically. Unfortunatelly, official documentation only provides the way how to do it through UI and does not provide corresponding REST API. I do not want to read the whole GeoServer documentation from ground up, just to learn at the very end, that it's imposible to do this task through API. So, I'm asking this question in advance with a hope, that someone can justify, that we can do all this - create a tiled layers - just using GeoServer REST API and without touching its UI. And to be more concrete, I guess, I need these three operations, addressed in the documentation:

creation of the coverage store by pointing at the directory containing the pyramid
choosing a coverage, which will be named after the pyramid root directory
publishing the layer, and then seting the layer parameter USE_JAI_IMAGEREAD to false to get better scalability


texts: Seeking shapefile of Indian electricity grid?

texts:
Does anybody know where I can find a shapefile of the major components of the Indian electricity grid?

Ideally, I'd get transmission/other high-voltage lines, substations, generating units, and transformers. I know there are power maps out there (like this one), but I want to match these transmission lines up with geographical data, so more precise locations would be helpful.


texts: Create continuous line network using PostgreSQL

texts:
I have shapefile which doesn’t have continuous network. I want to connect it to the nearest line and create continuous line from it. the postgis table has three columns- gid,objectid and geom.The geom is in multi line.now i want to connect this line to the nearest line and make it one line .
i try this query but its not connecting to the nearest line.

select st_makeline(st_linemerge(geom::geometry)) from footpath

.


texts: Difference between two layers in PostGIS

texts:
What is the right way to calculate the difference between two layers? I tried to use the following approach:

SELECT ST_Difference(river.geom, lakes.geom)
FROM river LEFT JOIN lakes ON ST_Intersects(river.geom, lakes.geom)

But in the output, I lose the geometries of the river layer which do not intersect any geometries in lakes. It looks like left join does not work as expected.

Currently, I'm using another approach, but I'm not sure this is correct:

SELECT ST_Difference(river.geom, lakes.geom)
FROM river JOIN lakes ON ST_Intersects(river.geom, lakes.geom)
UNION
SELECT river.geom
FROM river JOIN lakes ON NOT ST_Intersects(river.geom, lakes.geom)


texts:
I have a PostgreSQL database wich containes a table with type point an 4  tables without geometry.

I created foreign keys to garantee the relations between those tables. I migrated my db via postgis to QGIS. I made the the relations in QGIS once again.

After that I created a form where I wanted to add a relation in a tab.

 when I edit the feature where i added the relation I can see the fields in the tab.

Yet when I consult it from the attribute table when moving to form view I can't find it and when I creat a package to QField I also don't see it

if you have any idea about why I can't see my tab tell me.


texts: How to get intersection point from two lines (two coordinated points and two bearings)

texts:
I have tried for ages to find the algorithm for getting the coordinate of two intersecting lines from two given coordinated (latitude and longitude) points.

As the picture below suggests, I have two points (A and B) with known coordinates. With these two coordinates, come two bearings, or directions. I am looking for a way to find coordinate of C, where the lines-of-sight intersect.

Note: I am going to build it into a Java / C# Application. If there exists a library that already contains the math, that would help too.


texts: How to offset polyline geometries (not symbology) using QGIS?

texts:
I have a polyline layer which has a road network digitized, up the center, in both directions however the polyline are sat "one on top of the other".

I can display the lines separately by "offsetting" through the symbology however what i need is to change the actual positions of the lines by an offset so that the lines are then separate ( not sitting on top of each other)for exporting as a KML.


texts: How to run sextante algorithms outside of QGIS python console?

texts:
I'm kind of stuck trying to figure out the way to run sextante from a standalone python from OSGeo4W distribution. The reason I want to do this is that I got tired entering parameters in the dialog every time I want to test a model from Model Builder.

So here is the python script let's call it test.py

# as per http://qgis.org/pyqgis-cookbook/intro.html#using-pyqgis-in-custom-application
from qgis.core import *
# supply path to where is your qgis installed
QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
# load providers
QgsApplication.initQgis()

from sextante.core.Sextante import Sextante
Sextante.alglist()
Sextante.alghelp("saga:slopeaspectcurvature")

That I'm calling from my batch file

@echo off

set OSGEO4W_ROOT=C:\OSGeo4W
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%OSGEO4W_ROOT%\apps\qgis\python\plugins;%HOME%/.qgis/python/plugins
set PATH=%OSGEO4W_ROOT%\bin;%OSGEO4W_ROOT%\apps\qgis\bin;%OSGEO4W_ROOT%\apps\qgis\plugins

python test.py

The problem is that it says Algorithm not found whereas I get meaningful output from QGIS python console.

I feel like I'm missing to initialize something. But what?

Is there a better way to test a Model other than via entering tons of parameters using GUI?

UPDATE 7/2/2012

I'm looking for generic pythonic solution to test with "mine" algorithms. Aforementioned algorithm is just an example showing that something probably was not initialized.

UPDATE 7/27/2012

An alternative to Script Runner is to use IPython console to debug scripts. Other than that there does not seem to be a way to do simple unit testing with sextante with nothing else running:(

UPDATE 7/30/2012

As Victor Olaya suggests, I try to initialize Sextante like in the code below.

#!/usr/bin/env python

import sys
from PyQt4.QtGui import QApplication
from sextante.core.Sextante import Sextante

def main():
    """ main function or something """
    # as per http://qgis.org/pyqgis-cookbook/intro.html#using-pyqgis-in-custom-application
    from qgis.core import *
    import qgis.utils

    app = QApplication(sys.argv)
    # supply path to where is your qgis installed
    QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
    # load providers
    QgsApplication.initQgis()
    # how???
    # qgis.utils.iface = QgisInterface.instance()
    Sextante.initialize()
    run_script(qgis.utils.iface)

def run_script(iface):
    """ this shall be called from Script Runner"""
    Sextante.alglist()
    Sextante.alghelp("saga:slopeaspectcurvature")

if __name__=="__main__":
    main()

However I get something like

Traceback (most recent call last):
  File "test.py", line 29, in
    main()
  File "test.py", line 20, in main
    Sextante.initialize()
  File "C:\Documents and Settings\user\.qgis\python\plugins\sextante\core\Sextante.py", line 94, in initialize
    Sextante.addProvider(GrassAlgorithmProvider())
  File "C:\Documents and Settings\user\.qgis\python\plugins\sextante\grass\GrassAlgorithmProvider.py", lin
e 17, in __init__
    self.actions.append(DefineGrassRegionAction())
  File "C:\Documents and Settings\user\.qgis\python\plugins\sextante\grass\DefineGrassRegionAction.py", li
ne 16, in __init__
    canvas = QGisLayers.iface.mapCanvas()
AttributeError: 'NoneType' object has no attribute 'mapCanvas'

Well... it all becomes a mailing list discussion alike. Perhaps it worth moving to qgis-user or qgis-developer instead of SE.


texts: Can I create a custom Shapefile for SAS with QGIS?

texts:
Okay - I am in the process of configuring GIS within our organisation and am waiting for QGIS to be installed.

My aim is to use QGIS to compile a custom shapefile of the business regions used by our organisation that is 'snapped' to UK Postcode District (eg AB11) level.

In my mind (and from what i've researched) i think i will be able to build a UK layer on QGIS as i have received CODEPOINT with POLYGONS from Ordnance Survey.

I also know the postcode districts that make up each business region. So i could ask QGIS to 'dissolve' the postocodes into aggregated business regions based on this lookup table?

Then, my plan is to output a shapefile of the business regions to import into SAS.

As i've not got QGIS yet - am i too far from reality expecting the above?

Also is anyone up on using Shapefiles in SAS? If so, i have another question for you....

Thanks in advance.


texts: How to get QGIS to play with pandas (Windows)?

texts:
has anybody been successful getting pandas to work with QGIS's Python installation, on Windows? I've used the OSGeo4W installer to install QGIS2.2, which comes with Python 2.7 and Numpy 1.5.1. Any reasonably new version of pandas requires Numpy 1.6.1 or newer. So that seems like my main problem. I've tried:
- Building an older version of pandas (e.g. 0.7.3 and 0.4.1) using easy_install, but that gives me a very long list of errors, at least the last of which are related to distutils\msvc9compiler.py
- Using binaries of old pandas versions, but they are all built for Numpy 1.6.1
- Upgrading Numpy, which also gives me an error message. In case anybody is interested, the message ends with "File “C:\OSGeo4W\apps\Python27\lib\distutils\msvc9compiler.py”, line 299, in query_vcvarsall” newline “raise ValueERror(str(list(result.keys())))” newline “ValueError: [u’path’]"

I'm wondering if anybody else has ever been successful in getting pandas to work with a QGIS install on Windows??

It would be so easy if QGIS would ship with, say, a version of Numpy that's less than three years old. Does anybody know if such an upgrade is in the works for upcoming QGIS versions?


texts: How to update polygon field with point count for multiple rows? (Temporal point-in-polygon)

texts:
I have a shapefile with several regions. It contains basic info like ID, Code, Name1, Name2, etc.

It is overlaid with a shapefile of points (nearly 200k). Each point represents an event with a type and a day (YYYYMMDD). The table for the point layer looks like:

ID  DAY        STAT1  STAT2   STAT3 ....
01  20010901   23.45  53.12   12.34 ....
02  20010901   43.35  13.56   31.84 ....
03  20010902   35.01  21.61   11.47 ....

What I would like to do is:

Iterate through the point data for each day (one day may have one or several rows).
For each day, count the points in each polygon region.
Update the polygon table with a new column for that day with the point count as its value.

Ideally then, the final result would be a polygon feature with a table that resembles:

ID Code Name1 Name2 01 02 03 04 ... 198721
01 abc  ABC   A01   14  2 86 12 ...     27
02 xyz  XYZ   XZ7   66 47 23 88 ...     19

I'm using QGIS and believe this can be done w/ Python but I am not sure a) exactly how b) if another method might be more ideal.

Any help, pointers, or suggestions are greatly appreciated!


texts: How to get public access to edit a layer in geoserver

texts:
I want to develop a GIS system with OpenLayers and GeoServer which is installed with OpenGeo Suite 4.1 and I need to edit a layer with no GeoServer authentication. This means that guest members should have permission to edit a layer.

How can I do this?


texts: I was searching the web for benchmarks (but couldn't find any) that would indicate hardware performance (probably cpu in particular) in different GIS software in time consuming tasks (like hillshade generation, raster reproject, topology buildup), not for the purpose of comparing different GIS software but to compare performance for different cpus, multithreading impact, ssd, etc. If i started with that, I suppose i'd need a standarized set of data, publicly available so anyone could try and post results. But does anyone know if such comparisons already exist anywhere?


texts: "KeyError: location" for AOI Tile download in Planet Labs

texts:
I am using the download.py program to recreate the planet labs example to download all data in an AOI for specific dates. I am curious to understand where exactly do I insert my API Key
for now I added an additional line
os.environ['PLANET_API_KEY']="Actual API Key"
and
I am able to replicate all the steps except when I use
python download.py --download /tmp redding.json REOrthoTile visual

It gives me

**File "download.py", line 175 in module
    args.asset,args.overwrite)
File "download.py", line 124, in process download
   download url=result.json()[asset_type]['location']
KeyError: location"**

so I am unable to download any of the activated scenes.

Also is there a way to know if I am downloading analytic versus the RGB data.


texts:
I've successfully installed ArcSDE and ArcGIS for Server for 10.1. I'm running SQL Server 2012 on a Windows 2012 Server x64 machine.

I'm able to direct connect to the SDE geodatabase and version the data. So I assume this means that ArcSDE was installed and working correctly?

But I'm a little confused about the ArcSDE services and what exactly they are used for. I feel like I may not be using ArcSDE for what it was made for without using these services. Is this correct?


texts: Counting and print number of points from a CSV in a series polygons / buffered points in GEE?

texts:
I'd like to compare the number of calls from gps coordinates to the population around those coordinates in Google Earth Engine. My plan was to count the number of incidents (points from a table: Ola_Orders_13_May_2020) that fall within a series of polygons or buffered points and compare this to the population under the same buffered points and to export the results in a table.

I've summed the population within the buffered points from 'CIESIN/GPWv4/unwpp-adjusted-population-count', and printed as a table. I'd like to compare the number of occurrences (points) from the imported CSV table (Ola_Orders_13_May_2020) with the population in that area and print in a similar table.

My script below (missing the imported table):

//Map the Reports
var Incidents = table2.map(function (feature) {
  return feature.buffer(50, 1)
})
Map.addLayer(Incidents.draw({color: '000000', strokeWidth: 0}), {opacity: 1}, 'Reports')

var suburbs = ee.FeatureCollection([
ee.Feature(ee.Geometry.Point([46.688283, 24.598619]), {name: '23593'}),
ee.Feature(ee.Geometry.Point([43.921886, 26.391658]), {name: '19358'}),
ee.Feature(ee.Geometry.Point([46.732708, 24.691544]), {name: '15856'}),
ee.Feature(ee.Geometry.Point([39.825922, 21.4192]), {name: '10210'}),
ee.Feature(ee.Geometry.Point([55.3724, 25.277]), {name: '18934'}),

]).map(function(suburb) { return suburb.buffer(1000); });
print(ui.Chart.image.seriesByRegion({
    imageCollection: ee.ImageCollection('CIESIN/GPWv4/unwpp-adjusted-population-count'),
    regions: suburbs,
    reducer: ee.Reducer.sum(),
    seriesProperty: 'name',
  })
  .setChartType('Table')
  .setOptions({
    title: 'Estimated Population by location',
    vAxis: {title: 'Population in 1 Kilometer radius'},
  }));

https://code.earthengine.google.com/8e8439dbe7722237fb721f8c6cd8b99e


texts: ArcGIS Visibility  / 0 when observer point in raster cell center

texts:
As I discovered, the ArcGIS Visibility tool seems to produces a vertical  line from the observer point northwards in the above ground level raster and a 0 value in the visibility Raster, if the observer point is situated directly in the center of a raster cell (in my case derived from it).

Can anyone confirm this, or maybe does someone have a solution?

Moving the point and using a different raster resolution works, but is not what I intended.
I tried using Point and Point z; using SPOT and OFFSETA; setting said z values <, > and equal to the raster value.

I'm using ArcMap 10.4.1

edit: cell Values in the example are:

81,510223 | 80,945107 | 81,222015
82,140327 | 81,309288 | 80,910629
81,483345 | 80,921204 | 80,979355

Observer height is derived from the Raster. I also tried using fictious numbers greater and smaller than the raster value.

Interestingly this occurs only for single selected points or for more, if they are in vertical alignment, but not in horizontal alignment.

edit2:
After trying to find out, if this would affect my model, I discovered, that by using the extent environment, the fault is averted.
Still, I'm wondering what might cause this.


texts: Copying attributed relationship with ArcPy?

texts:
Is there a way using arcpy to copy an attributed relationship and only the attributed relationship ?

I can recreate a non attributed relationship with:

arcpy.CreateRelationshipClass_management(arg1,arg2,...)

From the doc this function support the argument "attributed" that can take the value "NONE" or "ATTRIBUTED", ok great, but now how to populate this relationship with the two primary key of the two feature class involved ?

I don't want to use arcpy.CopyFeatures_management because if we use this function to copy a relationship the function also copy the feature class involved.


texts: Save spatial reference to application settings?

texts:
Having created a tool as a macro for several geoprocessing tasks, I'm saving user input to application settings so the user doesn't have to enter all required data on each use.  I'm stuck on saving the spatial reference.

ESRI code shows using the FactoryCode and an enumeration to create either a projected or a geometry spatial reference.  What I'm trying to figure out is how to get the enumeration, and how to identify whether the factorycode/enumeration points to a projected or geometry spatial reference.

I will greatly appreciate any answers/directions provided.


texts: Mapping multiple routes from the same csv data file

texts:
I'm plotting a set of marches (each span 2-4 locations), and I've got all the latitude/longitude for each march within the same csv file. For marches with more than 2 locations, I've added an order ID field in the data, I've also added an ID field in the data that groups each of the marches.
I saw that Points2One and Points to Path would allow me to draw out multiple routes from the same set of data, but then I realised that it's not supported for the version of QGIS I'm using - 3.0.
I'm relatively new to QGIS, and I'm hoping someone might have suggestions as to how I can go about drawing these routes/what plug-in I could use so that I don't have to go in and manually plot each route.


texts: Clip raster not working: no target filename specified QGIS

texts:
I know this topic has been covered before but I didn't find an answer to my problem so far...
I've been trying to clip a raster unsuccessfully with canvas extent or with mask layer. It ever works and I always get the same message: "no target filename specified" or "no target data specified". With the specific message

*Loading resulting layers The following layers were not correctly generated.C:/Users/LongNhi
Computer/Documents/1.tifYou can check the 'Log Messages
Panel' in QGIS main window to find more information about the
execution of the algorithm. QGIS version: 3.12.0-București FAILURE: No
target filename specified.
Usage: gdalwarp [--help-general] [--formats] [-s_srs srs_def]
[-t_srs srs_def] [-to "NAME=VALUE"]* [-novshiftgrid] [-order n |
-tps | -rpc | -geoloc] [-et err_threshold] [-refine_gcps tolerance [minimum_gcps]] [-te xmin ymin xmax ymax] [-tr xres yres] [-tap]
[-ts width height] [-ovr level|AUTO|AUTO-n|NONE] [-wo
"NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16] [-srcnodata
"value [value...]"] [-dstnodata "value [value...]"] -dstalpha [-r
resampling_method] [-wm memory_in_mb] [-multi] [-q] [-cutline
datasource] [-cl layer] [-cwhere expression] [-csql statement]
[-cblend dist_in_pixels] [-crop_to_cutline] [-of format] [-co
"NAME=VALUE"]* [-overwrite] [-nomd] [-cvmd meta_conflict_value]
[-setci] [-oo NAME=VALUE] [-doo NAME=VALUE] srcfile*
dstfile*

My manip is the following, I'm a beginner in QGIS and I've been trying to solve it but nothing works. Anyone has an idea about what could I do?


texts: Division by zero error in Erdas imagine 2014

texts:
I'm running a Model maker to divide a couple of rasters. The function is as follows:

($n1_band1 / $n2_band2)* 100

There is no zero in both images (both of them are integers). However, I'm still getting the division by zero error.

Anyone know how I can solve this problem?


texts: Is there a way of creating a flowchart in QGIS or an alternative?

texts:
Currently I have thousands of animal trails and I want to analyse the network to see which waterholes trails lead to or if they are fragmented by roads. I tried exporting the CSV file to gephi another software that analyses networks but it doesnt create a flow chart. So I'm wondering if I could do a flowchart on QGIS or an alternative. Below is a partial of my project

And was aiming to get something like this...


texts: "libintl-9.dll is missing" when trying to run PostGIS 2.0 Shapefile and DBF Loader Exporter

texts:
When I try to run PostGIS 2.0 Shapefile and DBF Loader Exporter I get the following error:

  The program can't start because libintl-9.dll is missing from your
  computer. Try reinstalling the program to fix this problem.

I this previous question which asked a similar question about libintl-8. Following Devdatta Tengshe's answer I checked "C:\Program Files\PostgreSQL\11\bin" and the file is there:

then I checked pgadmin4's binary filepath and it's pointed at the same directory:

Where do i got from here?


texts: What is the purpose of featureNS in Geoserver?

texts: I wonder why geoserver have featureNS field, if it can not be retrieved by any publicly made request. AFAIK generally featureNS is company URL plus project name.


texts: Create road speed limit segments with linear referencing

texts:
I have a route-feature class and a point feature class with speed limit signs. In ArcMap I am capable to reference the signs on the route with help of a reference field and using the locate features along a route tool . So now I have the signs referencend on the routes and these referenced signs have a measure-field.

Now I am trying to create a line or route feature class with speed limit information on it but I can not find out how? The goal is that I have route segments that say: on route 1, from 0m to 1000m speed limit is 50kmh, from 1000m to 1337m speed limit is 80kmh and so on...
Any help is appreciated!


texts: Problem with Distance Matrix in Qgis 2.83

texts:
I want to calculate the distance of archaeological sites to the nearest river. Following this question
Calculating minimum distance between points and polygons in QGIS?
I extracted the nodes of the rivers layer and assigned a unique ID to each resulting point. I have a layer of almost 90,000 points/nodes. My archaeological sites layer has only 42 sites. Both layers are in EPSG 32615 (UTM 15N WGS84). I performed the Distance Matrix with the following settings:

Input Layer: Sites

Unique ID input layer: Site ID

Target Layer: nodes

Target unique ID: nodes id

Use only the nearest (k) target points" as 1

Now the problem, the resulting distance matrix states that all the sites have the same neighbor and are at the same distance of it:

This is not the only weird thing. The resulting distance should be in meters, since my layers are in EPSG 32615, but the resulting distance is almost two million meters!!. With Qgis manual measuring tool the distance between the “master node” to the nearest sites is almost 25,000 meters, and to the furthest one  about 56,000 meters.

I wonder if someone could give me a hint on what I am doing wrong.

Ps. I also tried MMqgis hub distance tool, but some of the resulting lines do not lie in a line node,

and the calculated distance is the same for every line


texts: Using pgr_drivingDistance for isochrones?

texts:
I have tried to calculate Isochrones in 500m walking distance along paths around planned subway stations with pgrouting.

For the setup I used an OSM-road.shp clipped to the borders of Vienna and a point shp. with the 14 stations which I want to analyse.

I worked my way through several tutorial since I am new to pgrouting for setting up a street-network and finding the closest nodes from the stations.

So far so good…but when I to calculate the isochrones with pgr_DrivingDistance the output shows very few results.

Either I set up the network wrong or I have been using the DrivingDistance function wrong. It seems like the function is not able to use every path.

How do I improve the query?

PostgreSQL 9.5.11
Pgrouting 2.4.1

CREATE OR REPLACE VIEW road_ext AS
  SELECT *, ST_StartPoint(geom), ST_EndPoint(geom)
   FROM road;

CREATE TABLE node AS
 SELECT row_number() OVER (ORDER by foo.p)::integer AS id,
    foo.p AS geom
FROM (
    SELECT DISTINCT road_ext.ST_StartPoint AS p FROM road_ext
    UNION
    SELECT DISTINCT road_ext.ST_EndPoint AS p FROM road_ext
   ) foo
  GROUP BY foo.p;

I had to change the geometry type since OSM provided multistring

ALTER TABLE road
 ALTER COLUMN geom TYPE geometry(LineString, 31256)
 USING ST_GeometryN(geom, 1);

Setup network

CREATE TABLE network AS
 SELECT a.*, b.id as start_id, c.id as end_id
  FROM road_ext AS a
    JOIN node AS b ON a.ST_Startpoint = b.geom
    JOIN node AS c ON a.ST_EndPoint = c.geom;

ALTER TABLE network add column len_km double precision

ALTER TABLE network add column time_walking double precision

UPDATE network set len_km = (ST_Length(geom))/1000

UPDATE network set time_walking  = ((60/5::double precision)*len_km)

Setup the point shp

ALTER TABLE ubahnpointsmgi
  ADD Column nearest_node integer

--Node in clostest distance is enough accuracy
CREATE TABLE temp AS
 SELECT
   DISTINCT ON (a.gid) a.gid, b.id AS start_id,
    ST_Distance(a.geom, b.geom) AS dist
    FROM ubahnpointsmgi a, node b
    ORDER BY a.gid, ST_Distance(a.geom, b.geom), b.id

UPDATE ubahnpointsmgi
 SET nearest_node =
  (SELECT start_id
    FROM temp
    WHERE temp.gid = ubahnpointsmgi.gid)

Calculating Line geometry

CREATE TABLE ISO_500 AS
 SELECT seq, cost, agg_cost, node, edge, l.geom FROM pgr_drivingDistance(
        'SELECT gid as id, start_id as source, end_id as target, len_km as
         cost FROM network',
        41217,0.5, false) AS foo JOIN network l ON l.gid = foo.edge

Calculating Point geometry

CREATE TABLE ISO_500 AS
  SELECT seq, cost, node, edge, l.geom FROM pgr_drivingDistance(
        'SELECT gid as id, start_id as source, end_id as target, len_km as
         cost FROM network',
        41217,0.5, false) AS foo LEFT JOIN node l ON l.id = foo.node

Output Table:
Output QGIS:


texts: ArcMap Mosaic Dataset function chain result to NetCDF

texts: I've been working on a mosaic dataset and added some changes by using the function chain options. This works fine in ArcMap and I can create a time varying result out of this. Since I want a 3D result as well using ArcScene and ArcScene can't handle mosaic datasets, I tried to convert the Mosaic Dataset to NetCDF. Unfortunately it seems that ArcMap doesn't convert the final result from the function chain. How can I get this to work, so I can get my function chain result visualized in ArcScene?


texts:
I have the same question as described in this question.

I have numerous milepost addresses (like Milepost 58.3 , US 30,NJ)

Is there a way to do this in PostGIS?


texts: I've been trying to fill those lidar scan lines in my dem with no avail. Is there any way to interpolate those "gaps" to get a clean dem? I used las2dem to create the dem from .laz file


texts: What makes a MXD file size larger and how to decrease its size?

texts:
I know that MXD contains information about layers, symbology and configuration of default tools. I have a few that are exceeding 50 MB and are starting to take a while to load up.

What are some of the main things that will increase a MXD's file size? Is there something specifically that dramatically increases the size?

Some of my ideas, that I have absolutely no proof of are:

Complex layer symbology
Multiple label classes per layer
Definition queries
Pulling layers from sources will long file paths (C:/MyFolder/Data/Data2/SubContainer/YouGetThePoint.shp)
Scale Ranges for layers
Arranging layers into groups

Do any of these ideas have any validity?


texts: Python QGIS, attribute values are saved as NULL using actionAddFeature().trigger()

texts:
I am new using python in QGIS. I was trying to create a simple plugin in QGIS to create a new shapefile with a single attribute 'id', start the edition mode, and add new features manually on the canvas.
shp_output_path = 'C:/test/P8_100.shp'

layerFields = QgsFields()

writer = QgsVectorFileWriter(shp_output_path, 'UTF-8', layerFields, QgsWkbTypes.Point, QgsCoordinateReferenceSystem('EPSG:25830'), 'ESRI Shapefile')

editlayer = QgsVectorLayer(shp_output_path, 'P8_100', "ogr")

layer_provider = editlayer.dataProvider()
layer_provider.addAttributes([QgsField('id', QVariant.Int)])
editlayer.updateFields()
del(editlayer)

#I reload the layer, otherwise I could not add new features
editlayer = QgsVectorLayer(shp_output_path, 'P8_100', "ogr")
QgsProject.instance().addMapLayer(editlayer)
editlayer.startEditing()
iface.actionAddFeature().trigger()

So, when I click on the canvas the usual window to insert the id value pops up. I insert a value. After inserting all the values, and before saving changes, I open the attribute table, and I see that the values inserted appear in the id attribute of the attribute table. However, if I save the changes, these values automatically change to NULL.
This does not happen when I launch this code from the python console of QGIS. It happens when the code is launched as a plugin.
I would like that the id values inserted in the pop-up window are inserted in the attribute table after saving changes manually or using commitChanges() command.


texts: Error when trying to populate a Dictionary with arcpy.da.SearchCursor using file paths and field name lists

texts:
I have a dbf table with two fields of interest. Field one is named AP_ZIP and rows contain field names for a zip code field within shapefiles. In addition, Field two is named APFI_PATH and contains the file path to a shapefile. I am using these two fields as parameters for arcpy.da.SearchCursor

In the code below, I created a list of file paths for shapefiles called filepathlist and a list of corresponding field names within the shapefile called zipfieldlist.

I'm trying to feed the file path from filepathlist and corresponding field name from zipfieldlist into the arcpy.da.SearchCursor and populate a dictionary. I want the dictionary to have file paths as keys and row items from fields as values.

import arcpy

arcpy.env.workspace = r"C:\AP_PAR"

dbf = "APExists_Test.dbf"
appath = 'APFI_PATH'
zipfield = 'AP_ZIP'

filepathlist =[row[0] for row in arcpy.da.SearchCursor(dbf,'APFI_PATH')]
print(filepathlist)

zipfieldlist =[row[0] for row in arcpy.da.SearchCursor(dbf,'AP_ZIP')]
print(zipfieldlist)

ZipDict = {}

for fp in filepathlist:
    for fields in zipfieldlist:
        with arcpy.da.SearchCursor(fp,[fields]) as curs:
            for row in curs:
                keyvalue = fp
                ZipDict[keyvalue] = [row[0]]
                #print(row[0])

print(ZipDict)

When I run the code I get the following error.

However, if I uncomment print(row[0]), the code runs and starts to print out zip codes. In addition, the dictionary ZipDict gets populated with one and only one key:value pair {'C:\\AP_PAR\\AL\\tuscaloosa\\E911_GDB_SDE_address.shp': [('35404',)]} len = 1

My desired output is to have a dictionary with keys as APFI_PATH and values as AP_ZIP. I want the values to be a list of strings. Ideally I want the values to be integers and not strings, but one step at a time.

Desired Output Example: Key = APFI_PATH and Value = AP_ZIP (list of strings)

{'C:\AP_PAR\AL\tuscaloosa\E911_GDB_SDE_address.shp': ['35404','35401','35402', etc....]}

Any ideas as to why I am getting this error and how to fix the code?


texts: Unable to install pgrouting for postgres 9.3.5/postgis 2.1.7

texts:
I am running postgres 9.3.5 and postgis 2.1.7 on Ubuntu Server 14.10. I am trying to install pgrouting, but having no luck.

1) I tried the following based in the official workshop documentation:

sudo apt-add-repository -y ppa:ubuntugis/ppa
sudo apt-add-repository -y ppa:georepublic/pgrouting
sudo apt-get update
sudo apt-get install postgresql-9.3-pgrouting

Then I get this error:

E: Unable to locate package postgresql-9.3-pgrouting
E: Couldn't find any package by regex 'postgresql-9.3-pgrouting'

2) I tried using the UbuntuGIS repo as well, but had no luck.

Any idea how to install pgrouting on my setup? I suspect that this is a Utopic Unicorn issue, but not sure.


texts: Most efficient way to rasterize a planet

texts:
I'm working on a Webservice where GIS supports my services as a Crowd Intelligence Aggregator to help civilian (like civilized) aid to support provisioners in a case of emergency, e.g. pestilence, military/civil conflict or other natural planet-bound catastrophies like earthquakes.

Example:
I imagine, that a convoy tries to support folks with vaccines. So it has to decide, whether it drives one road or a parched stream bed. My planned open source service will tell via visible heatmap, whether it may pass easily to its destination via stream bed, or it's likely to be attacked by brigands (where it happened the past weeks) on the road though noticed and announced to the webservice.

To rasterize the planet I'm thinking of how to store local gathered information on a map/datawarehouse as a synonym for fast predictive information analysis. The tiles shall be a tetragon or pentagon. The goal is with regard to calculate a - heatmap. Maybe it's harder to tesseract a pentagon than an tetragon... Maybe you folks have experience or knowlegde to help me? I'm not sure what's "better"...


texts: Grass GIS r.to.vect: Prevent neighboring cells with equal values from being unified

texts:
I've got different satellite image bands as GRASS rasters, which I want to convert to vector dataset. As I want to compare the IDs of the cell polygons later on, the neighboring cells with coincidentally equal values should NOT be unified to one polygon like this:

Is there a way to prevent this case? I can't use the '-v' flag to use raster values as categories since I just need the unique IDs later on.


texts: Cannot find python executable being used by QGIS

texts:
I am trying to develop a plugin for QGIS but I am having trouble setting up my IDE with the required python executable so I can get a console working for logging output while developing it.  When I deliberately enter some messed up text in the script, and then try to run the plugin in QGIS, I can see an error in the stacktrace that I would assume indicates where the executable is:

Traceback (most recent call last):
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 333, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/isaac/.qgis2/python/plugins/MakeOSMRoutableNetwork/__init__.py", line 34, in classFactory
    from .make_osm_routable_network import MakeOSMRoutableNetwork
  File "/usr/lib/python2.7/dist-packages/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/isaac/.qgis2/python/plugins/MakeOSMRoutableNetwork/make_osm_routable_network.py", line 214
    self.dlg.postgis_connection_listView.addItem(k)weqweqfasfsdf
                                                               ^
SyntaxError: invalid syntax

So it seems like my python exectuable is in /usr/lib/python2.7.  Alright, fair enough.  But when I navigate to this directory and list the contents, I don't see any python executable.  Shouldn't there be one in here? I am using Ubuntu 16.04.


texts: Language requesting GeoDjango technical documentation

texts:
I am working with a project team that hired a contractor to build a custom GeoDjango webmap application.  The application is mature now and the team would like to host it on our network "in-house".  I'm responsible for writing the contract deliverable to get the technical documentation (something that was overlooked during the building of the application, before I was hired here).

How much detail should I be looking for in the technical documentation?  Currently, the github site has some sparsely commented installation instructions and nothing else.  But, if a few months down the road I need to get under the hood to add functionality, how much documentation should I expect?  Aside from the full technical specifications, what are the key details that are must haves to facilitate a smooth transition?

I need to make sure this captures everything we'll need to know moving forward.  Right now I have a very generic:

Create basic training materials for new users in the form of a user manual as well as technical documentation providing full specifications of the technical stack…


texts: Stream network using multiple cell thresh holds inArcGIS

texts: I'm trying to create a stream network with drainagelines formed at varying thresh holds throughout a drainage basin. Is there any way of doing with Arc and still creating a stream network that is fully connected?


texts: I need to export DXFs with a certain Z value

texts: How do I do that? Right now it says it's including a Z dimension but I don't know where to input it.


texts:
I have a .csv file with population data according to specific locations (Points).
My problem is that all locations are given in WKB format like that:
0101000020E6100000B953EFEEEEB253C0F15BD4CCCC4C05C0

How can I visualize my data with QGIS?
I don't know why there is no simple solution. Because to visualize WKB from a database is no problem.
I hope you can help me.


texts: Is there any way to freeze the layer names in the TOC?

texts:
I would like to freeze the names of layers in TOC.
Is there any way to do this via QGIS config or PyQGIS?

I´m using QGIS Wien.


texts:
I was meandering so many places to find GIS data for Beijing, but it seems like I cannot get to any useful resources.  Though I did get at least road, water shapefiles and outline of Beijing city.

Specifically, what I need is buildings and parcel (lots) maps to investigate site analysis in the Beijing CBD area.


texts: Using Shapely & Python to remove overlapping polygons in one dataset

texts:
I am looking for a method using python to remove a large amount of overlapping polygons from one collection. These polygons are grouped together based on imagery collection and have been stored as WKB in a database. I am able to pull the WKBs out per collection. My current approach is to create a list of geometries and iterate using shapely.intersects per polygon to an empty list. If polygon in listA does not intersect with polygon in listB the polygon in listA is appended to listB and removed from listA. My thinking is as this iterates I will be able to move polygons with no overlaps to one list and keep the overlapping polygons in another. When I run this code I get the error: "AttributeError 'list' object has no attribute '_geom'".

polygon_list_A = [<shapely.geometry.polygon.Polygon object at 0x000005B3A>, <shapely.geometry.polygon.Polygon object at 0x000006E7A>, ...]

polygon_list_B = []

for polygon_A in polygon_list_A:
    if polygon_A.intersects(polygon_list_B) == True:
        polygon_list_B.append(polygon_A)
        polygon_list_A.remove(polygon_A)
    else:
        pass


texts: Changing the style of a polygon with a click event in a Shiny/Leaflet app

texts:
I'm attempting to create my first Shiny app using Leaflet. I'm currently trying to figure out how observers work together to update the map when something is clicked. Here is the example that I've been working from. In this example (code included in link), each marker changes styles and becomes yellow when clicked on, shown below:

Similarly, in my app, I want my polygons to change styles when they are clicked. While the example is very useful, I'm having problems translating the code from marker clicks to shape clicks.

Here is some example code with my shapefiles available for download here.

library(shiny)
library(rgdal)
library(leaflet)

ctry <- readOGR("D:/", layer = "ctry", stringsAsFactors = F)

shinyApp(

  ui = fluidPage(leafletOutput("map")),

  server = function(input, output, session) {

    output$map <- renderLeaflet({
      leaflet(ctry) %>%
        addTiles() %>%
        addPolygons(fillColor = "gray",
                    fillOpacity = 1,
                    weight = 2,
                    stroke = T,
                    color = "blue",
                    opacity = 1)
    })

    observeEvent(input$map_shape_click, {
      click <- input$map_shape_click
      proxy <- leafletProxy("map")
      if(click$id == "Selected"){
        proxy %>% removeShape(layerId = "Selected")
      } else {
        proxy %>% addPolygons(fillColor = "black", layerId = "Selected")
      }
    })

  })

This code pretty closely follows the linked example with an attempt to work with shapes instead of markers. When I run the code, however, I get the following error:

  Warning: Error in if: argument is of length zero Stack trace
  (innermost first):
      68: observeEventHandler [#22]
       4:
       3: do.call
       2: print.shiny.appobj
       1:  ERROR: [on_request_read] connection reset by peer

I'm not sure where to go from here. Again, I want to change the style of my polygons whenever they are clicked (in the example, I made the fill black). What am I doing wrong?

I've continued to play around with my code a bit more and I think I have narrowed the problem down to the if(click$id == "Selected")... statement. It was my understanding that the click event is what triggered click$id to be populated by the string "Selected." However, when I print click$id from the example markers, I get the location name. I am really unsure of how this "Selected" id works.


texts: How can I reclassify an integer raster to a float raster?

texts:
I have been using various arcGIS tools to try and do this. Unfortunately, the reclassify tool in arcmap can only output integers. A cool workaround that I found is using the raster calculator in arcmap, as suggested by Jeff Evans in this post: http://forums.arcgis.com/threads/15396-Reclass-from-Integer-to-Floating-point

*"Here is an example using a 4 value integer (21, 41, 42, 43) raster that results in a float with the values of 1.0, 2.0, 3.0, 4.0:
Float(Con("nlcd" == 21, 1.0, Con("nlcd" == 41, 2.0, Con("nlcd" == 42, 3.0, Con("nlcd" == 43, 4.0, 0.0)))))"*

This works like a charm and gives me the desired output...

HOWEVER, I cannot use this process in arcpy... Because I need to do many loops of this process, I need it to be its own python script and “The Raster Calculator tool is intended for use in the ArcGIS for Desktop application only as a GP tool dialog box or in ModelBuilder. It is not intended for use in scripting and is not available in the
ArcPy Spatial Analyst module”.

To be honest, I am starting to give up on using arcpy for this, does anyone out there know how to do this efficiently using arcpy? I also tried some R code but it took way too long (my grids are around 1500x2400 cells). I also tried a python reclass script my colleague he built scratch, and it also takes too much time.

Any ideas/suggestions?

Thanks for your time!


texts:
Does ogrinfo support a GROUP BY clauses? I could not find any good documentation of all the SQL syntax features it supports.

I tried the following:

ogrinfo . -geom=no -sql 'SELECT lsad FROM tl_2011_us_county GROUP BY lsad'

And got the following error:

ERROR 1: SQL: Failed to parse SELECT statement, extra input at BY token.
INFO: Open of `.'
      using driver `ESRI Shapefile' successful.

What I really want to do:

SELECT lsad, count(*)
FROM tl_2011_us_county
GROUP BY lsad

What I can do:

SELECT DISTINCT lsad
FROM tl_2011_us_county

And then:

SELECT count(*)
FROM tl_2011_us_county
WHERE lsad = 'foo'

... but that is annoying.


texts: Rasterio: Issue with `writing a dirty block from GDALRasterBand::IRasterIO` even with `LibTiff` > 4.0 on Windows and Mac, but not Linux

texts:
I was hoping someone might have some insights. I am using python rasterio to compress some large GeoTIFF imagery files. By large I mean 2-3 GB compressed, or on the order of 60,000 x 60,000 pixels each. I am hitting an error with BIGTIFF that should have been resolved since I am using Libtiff greater than 4.0, but the compression process is crashing on large files on Windows and Mac. So I followed the advice of previous posts on this topic.
I am using rasterio 1.1.7, and Gdal 3.1.2, and LibTiff 4.1.0. The code I am using for compression is below, and the error message is below as well. Note I am also using "windowed" writing to help avoid loading all of the data into memory--so I should not be hitting memory constraints that way.
dat = 'sample_image_big.tif'
target_file = 'sample_image_big_compressed_jpeg.tif'

with rasterio.Env():

    profile.update(
        compress="JPEG",
        BIGTIFF = "YES")

    with rasterio.open(target_file, 'w', **profile) as dst:
        for ji, window in dat.block_windows(1):
            dst.write(dat.read(window=window), window=window)

And here is the corresponding error message:
Traceback (most recent call last):
  File "rasterio/_io.pyx", line 1369, in rasterio._io.DatasetWriterBase.write
  File "rasterio/shim_rasterioex.pxi", line 133, in rasterio._shim.io_multi_band
  File "rasterio/_err.pyx", line 182, in rasterio._err.exc_wrap_int
rasterio._err.CPLE_AppDefinedError: An error occurred while writing a
dirty block from GDALRasterBand::IRasterIODuring handling of the above
exception, another exception occurred:Traceback (most recent call
last):
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/bin/nyuki", line
10, in <module>
    sys.exit(nyuki())
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/click/core.py",
line 829, in __call__
    return self.main(*args, **kwargs)
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/click/core.py",
line 782, in main
    rv = self.invoke(ctx)
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/click/core.py",
line 1259, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/click/core.py",
line 1259, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/click/core.py",
line 1066, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/click/core.py",
line 610, in invoke
    return callback(*args, **kwargs)
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/nyuki/nyuki.py",
line 160, in compress
    compressor(sourcetiff, target_compression, yes)
  File "/Users/user/mdconda/anaconda3/envs/nyuki-env/lib/python3.7/site-packages/nyuki/geotiff_compressor.py",
line 52, in compressor
    dst.write(dat.read(window=window), window=window)
  File "rasterio/_io.pyx", line 1371, in rasterio._io.DatasetWriterBase.write
rasterio.errors.RasterioIOError: Read or write failed. An error
occurred while writing a dirty block from GDALRasterBand::IRasterIO

If anyone wants the actual test image, they can pull it at this link. I am not sure what the problem is? I am using a "windowed" read and write to conserve memory. The hard disk has 100 GB of space, so no problem there. So I am not sure what this error is indicating or how to fix it. Any help is appreciated.


texts: Error with z values in ArcGIS ModelBuilder

texts:
Although the input data has z values the tool "Adjust 3D z values" in ModelBuilder gives

  ERROR 000245: The input feature class must have z values

A multipoint feature gets intersected with a polygon feature and then the intersected multipoints get elevated with "Adjust 3D z value".
I have already set the properties of every feature in the model to "z-values: Enabled" or "z-values: Same as Input" but this doesn't help. As the model is very big it is not possible to post a readable screenshot.

I'm using ArcGIS 10.2 for Desktop Advanced with 3D-analyst.

Any ideas what the problem could be?


texts: Reading bathymetry for coordinate with Python?

texts:
I need the water depth at specific coordinates and my geoscientist are offering me a dataset in either TIFF or XYZ/txt format with data points each 2km.

I will typically be extracting about 1000000 points from the dataset in a single run of my application, and the dataset covers the entire Norwegian Continental Shelf.

I have never dealt with spatial data before, so how do I extract the depth a certain x,y coordinates?

The data would be used within an application written in Python, but of course a backend structure that parse the dataset could be in any other language as long as it provides an interface I can use in Python.


texts: Generating OSM map without HTML rendering but PDF rendering with iText?

texts:
Is there a way to create an OSM map by a java bean, without displaying it in a web page, but instead stamping it inside a PDF which is generated with iText 5?

In my case, I want to generate kind of a serial letter as PDF, each copy containing a map, that displays a personal routing for the recipient.

So there is no need to show a map in a jsp or JSF page, only to create it and add it to the PDF.


texts: Use line style to adjust rubberband in QGIS?

texts:
I am digitizing lines and have to keep certain distances to other lines while digitizing.
Normally I would just use a buffer to generate these prohibited areas and snap to these polygons while digitizing. But I have only access to wms-layers so I cannot use snapping.
Therefore I would like to adjust the rubberband so I could always see a 5,10 and 20m distance-(rubberband)line while digitizing.

Is it possible to adjust more than just the color for the rubberband so I could always see these distances while digitizing?
Here is a example screenshot to show how it could look like (dotted line=digitized line, other three lines = distance lines).


texts: Update an empty column with a geometry selection (ST_Centroid) from the same table (or another one)

texts:
I'm admitting that PostGIS is really giving me headache and it's not my strong point.
I have a table my_table with columns: id_point , point (point n°, contains numbers from 1 to 5), geometry (geometry of each point), road (road number), geom_moyenne.
Geometry has a point geometry for each point, geom_moyenne is empty.
I have made a selection to get an average geometry for each road (each road has 5 points with different geometries and I want to have a kind of a middle position point for each road):
SELECT ST_Centroid(ST_Union(geometry)) AS geom_moyenne , road
FROM my_table
GROUP BY road;

Now I want to put the results of my selection in my empty geom_moyenne column and I failed miserably. Neither of these two is working. What am I doing wrong?
INSERT INTO my_table (geom_moyenne)
SELECT ST_Centroid(ST_ Collect (geometry))AS geom_moyenne
FROM my_table
GROUP BY road ;

--
UPDATE my_table
SET geom_moyenne = ST_Centroid(ST_ Collect (geometry)) AS geom_moyenne
FROM my_table
GROUP BY road;

I also have a second table my_roads that has road (same road numbers as my_table)and geom_moyenne columns, and I'm trying to put the same selection in the geom_moyenne column, but can't do that either (I need it in my_roads table actually but just thought it would be easier to first fill it in the same table and then eventually copy to another).


texts: Determine if approximate count of features in polygons is less than 100_000

texts:
Using PostgreSQL 12.
I have a table of polygons, a table of points, and a table of lines. All tables have GIST indexes on the geometry column. There is a view that includes both the lines and points, view_points_and_lines.
There is a script that I run on these that is fairly intensive, and I do not want to run it if there are more than 100_000 total features (points and lines) in the targeted polygons (combined).
Thus, before I run the "intensive query", I would like to determine whether the total count of the features in the targeted polygons is < 100_000. Note that if a feature is included in multiple polygons, it should be counted multiple times.
However, running a query like
SELECT COUNT(tf.*) FROM view_points_and_lines tf JOIN polygons sf ON ST_Intersects(tf.geom, sf.geom)

can be very slow itself, especially if there many polygons and much more than 100_000 points.
I have also tried returning a Boolean like:
SELECT COUNT(tf.*) < 100000 FROM view_points_and_lines tf JOIN polygons sf ON ST_Intersects(tf.geom, sf.geom)

but that does not seem to speed up the query, and the query plan is the same.
Given that I really only need a Boolean value, and the count can be approximate, is there any way to speed up this query?


texts: .ecw file not imported in its projection in QGIS

texts:
I have a .ecw aerial image and one of the associated files is .prj, which i assume defines the Coordinate Reference System of the image. The .prj is as follows:

PROJCS["PNG94_PNGMG94_zone_55",GEOGCS["GCS_PAPUA_NEW_GUINEA_GEODETIC_DATUM_1994",DATUM["D_PAPUA_NEW_GUINEA_GEODETIC_DATUM_1994",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",147],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]

Which both here and here is described as EPSG:5551. Based on this I chose EPSG 5551 as the CRS of my QGIS project and made sure in settings that "When a new layer is created, or when a layer is loaded that has no CRS" QGIS "Use project CRS".

However, when i import my .ecw file, the project CRS is changed to EPSG 7855 (GDA2020 / MGA zone 55) and when i do right click on the layer and go to properties, i see that the layer's CRS is EPSG 7855. In other cases, another image's CRS changed to EPSG 28355.

I know that EPSG 5551, 7855 and 28355 are related, and all cover Zone 55 between Australia and Papua New Guinea. I want to understand why the CRS of my project is changing? Why the layer appears to have another CRS when i check its properties in QGIS? and if this may have a consequence in my project?


texts: Calculating euclidian distance in R between lines and points?

texts:
I have a large SpatialLinesDataFrame "o" (~500,000 roads) and ~3000 survey points "r", I'd like to calculate the distance to the nearest line for each survey point and extract relevant columns.

My solution at present is to split the SpatialLinesDataFrame into chunks and calculate distance using gDistance in the "rgeos" package. This will take over 18 hours though on my machine.

Is there a faster way to do this, perhaps by converting lines to points first?

chunk <- 10000
n<-dim(o@data)[1]
m<-rep(1:ceiling(n/chunk),each=chunk)[1:n]

m.list<-vector("list", length=dim(r)[1])
for (i in 1:dim(r)[1]){
df.list<-vector("list", length=length(unique(m)))
for (j in 1:length(unique(m))){
p<-o[m==j,]
df.list[[j]]<-data.frame(ID=i,Dist=min(gDistance(r[i,],p,byid=TRUE)),p@data[which.min(gDistance(r[i,],p,byid=TRUE)),])}
d<-data.frame(do.call("rbind", df.list),stringsAsFactors=F)
m.list[[i]]<-d[which.min(d$Dist),]
}


texts: Using UpdateCursor in the same loop as SearchCursor

texts:
I have 2 layers, one with land plots and one with hazard zones.
I want to iterate through the layer with land plots and for every plot I want to add 2 (or 3, depending on some variables) fields. Then I calculate the data for these fields, based on the attribute values of the overlapping hazard zone, and fill it in, in the new fields. I'm trying to do this by using the searchcursor to iterate through the plots and get their geometry (I need the geometry to decide if they overlap, touch, intersect etc. with the hazard zone, which affects the calculated data for the new fields).

Next to all this I try to create an updatecursor and loop through both cursors at the same time using zip(), to update the newly created empty fields immidiately after calculating their proper value.

I'm trying to do this in 1 for loop to save time by only iterating the plots once.
But the UpdateCursor doesn't work properly. The 2 new fields are Null for all plots except the last plot, that one has the proper values.

Any clues to what is happening here?

I have printed the variables which should be added and they are not Null.

Start of the loop:

        fields_og = arcpy.ListFields(outLayer)
        fields = ['SHAPE@']
        for field in fields_og:
            fields.append(field.name)
        with arcpy.da.SearchCursor(outLayer, fields) as features:

            if atype == 'SHORT':
                arcpy.AddField_management(outLayer, attrn, atype, 10, 0)
            elif atype == 'DOUBLE':
                arcpy.AddField_management(outLayer, attrn, atype, 10, 2)
            else:
                arcpy.AddField_management(outLayer, attrn, atype)

            attrn_s = attrn + "_s"
            minmax_s = attrn + "_mm"
            arcpy.AddField_management(outLayer, attrn_s, 'TEXT')
            out_fields = [attrn, attrn_s]
            minmax = False
            if len(a_minmaxfield) > 0:
                arcpy.AddField_management(outLayer, minmax_s, 'TEXT')
                minmax = True
                out_fields.append(minmax_s)

            with arcpy.da.UpdateCursor(outLayer, out_fields) as cursor:

                for feature, row in zip(features, cursor):

After some calculations etc. I want to update like this:

                    if minmax:
                        attrs.append(c)

                    if minmax:
                        row[-3] = val2add
                        row[-2] = b
                        row[-1] = c
                    else:
                        row[-2] = val2add
                        row[-1] = b

                    cursor.updateRow(row)


texts: QGIS failed with error: latitude or longitude exceeded limits

texts:
I have some shapefiles from vegetation types in QGIS and I wanted to put a Google map in the background, so I installed the OpenLayers plugin. But every time I tried to load the map, there's an error message:

"failed with error: latitude or longitude exceeded limits".

What to do then?


texts: What is the best way to sum the area from a number of polygons in .NET/SQL Server?

texts:
I have an Mapbasic/.NET application that read data from SQL Server 2008 R2. I have a table with polygons and need to sum up the area of a subset of the polygons. This should work dynamic so the user should be able to sort on different parameters and get the sum of the areas from their selection. My options are:

MapBasic
C# .NET
SQL SERVER

Mapbasic seems to be to slow for this, personally I would prefer to do this direct in SQL Server but have not found any way to sum up the area with T-SQL. Have anyone done this before or is there any reason to do this at any of the levels?


texts:
The following KML file contains paths with style and name tags: example.

<Folder>
          <name>Pedestrian</name>
          <Style id="style-Pedestrian">
            <LineStyle>
              <color>#ff808080</color>
              <width>3</width>
            </LineStyle>
          </Style>
          <Placemark id="">
            <styleUrl>#style-Pedestrian</styleUrl>
            <LineString>
              <altitudeMode>relativeToGround</altitudeMode>
              <coordinates>8.5950278,46.6329618,1 8.5951371,46.6329856,1 8.595266,46.6327811,1 8.5951712,46.6327501,1 8.5950278,46.6329618,1</coordinates>
            </LineString>
          </Placemark>
....

I tried to refer to the styles, but the TileMill does not recognize it. After I import the files and apply the style using "Carto styling language" the path just disappears.

Example:

#imported_id {
  [Name="Pedestrian"] {
    line-width:2;
    line-color:#000000;
  }
}

#imported_id {
  [Name="style-Pedestrian"] {
    line-width:2;
    line-color:#000000;
  }
}

#imported_id  #style-Pedestrian {
      line-width:2;
      line-color:#000000;

  }

Why I can't refer to the styles inside the TileMill? Is it possible or I am doing something wrong?


texts: 'Must have equal len keys and value when setting with an iterable' error message

texts:
I have a postgres spatial table that I want to convert to a geopandas GeoDataFrame. In order to do so I wrote the following piece of code:

#'pg_tab_clip' is the name of the postgres spatial table
#'b' is an integer
#'con' is the connection to postgres

    select = etl.fromdb(con, "SELECT ST_AsText(geom) FROM {0} WHERE buffer = {1}".format(pg_tab_clip, b))
    temp = select[1:len(select)]
    buffer_polygs = [shapely.wkt.loads(t[0]) for t in temp]

    buffer_gdf = gpd.GeoDataFrame()
    for i, bp in enumerate(buffer_polygs):
        buffer_gdf.loc[i, 'geometry'] = bp

When I run the script I get the following error:

Traceback (most recent call last):
  File "02_culvert_processing_v1.py", line 417, in <module>
    main(working_folder, alr_path, pg_tab_culv, N, margin, pg_con)
  File "02_culvert_processing_v1.py", line 390, in main
    loop_tiles(N, margin, con, cursor, alr_path, srid, pg_tab_buf, pg_tab_tile1, pg_tab_clip, buf)
  File "02_culvert_processing_v1.py", line 274, in loop_tiles
    process_raster(con, pg_tab_clip, tile2_tif, srid, i, j, margin, culv_count, buf)
  File "02_culvert_processing_v1.py", line 143, in process_raster
    buffer_gdf.loc[i, 'geometry'] = bp
  File "C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\indexing.py", line 190, in __setitem__
    self._setitem_with_indexer(indexer, value)
  File "C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\indexing.py", line 366, in _setitem_with_indexer
    self._setitem_with_indexer(new_indexer, value)
  File "C:\ProgramData\Anaconda3\lib\site-packages\pandas\core\indexing.py", line 611, in _setitem_with_indexer
    raise ValueError('Must have equal len keys and value '
ValueError: Must have equal len keys and value when setting with an iterable

I don't know what's causing this problem. The strangest of all is that the script doesn't actually fail every time I run it but only sometimes (???).
Also, is there perhaps a more elegant/pythonic way to do what I am trying to do? Is the iterable 'buffer_polygs' actually unnecessary?