texts: How do I set layer transparency in QGIS 2.0 with Python?

texts:
I have a raster layer loaded in QGIS and I'd like to be able to set its transparency via the Python console.  I'm new to the API and have not been able to find out how to do this.  Any help would be greatly appreciated!

Update:

Here's what I see from the console:

>>> l = qgis.utils.iface.activeLayer()
>>> l
<qgis.core.QgsRasterLayer object at 0x07144DF8>
>>> rt = l.rasterTransparency()
Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'QgsRasterLayer' object has no attribute 'rasterTransparency'

I'm just loading a single band raster in GeoTiff format with a palleted color scheme... nothing too fancy.  Am I getting the layer incorrectly?


texts: Midvatten Stratigraphy "wrong table" error

texts:
I am creating a dummy data set to do some geological modelling on QGIS using Midvatten, and subsequently Qgis2threejs plugin to create a 3D borehole visualisation.

Is there a missing "Stratigraphy" column between WQual and Piper?

And any idea why there's a "wrong table error?"


texts: ArcGIS Online app multiple comments on one feature

texts:
I am working on an ArcGIS application that would allow the public to comment on bus routes though a map based interface. The original thought was to use the group sort tool to easily navigate through the routes and use the edit tool to leave comments, but I believe that this would be a bit confusing for the public.

What I would like to accomplish is an application that when you clicked on the route the pop-up would give you the option to leave a comment either right in the initial pop-up or a link to an additional pop-up.

It seems like the application would be in between the edit application and the geoform application.

Has anyone built something like this before


texts: Database connection (SDE) failed after database was renamed

texts:
I created an enterprise geodatabase in SQL Server 2008 R2 with ArcCatalog (10.2). It was connected via an sde file in ArcCatalog then.

Then I had to rename the database from SQL Server Management Studio.

After I did renaming the sde file no longer worked. The message is:

  Failed to connect the specified server. Do you want to continue?
  Operation Failed

I am sure the connection properties are correct as I only rename the database, and the Database I selected to connect is the new name.

Any idea?


texts: Editing relationship class in ArcGIS for Desktop?

texts: I want to know if there's a way to edit the origin/destination class in Relationship class in arcgis after creating the Relationship class?


texts: Referencing a polygon layer in the field calculator of a point layer to make a virtual field

texts:
I have a polygon layer that shows the counties of England and a point layer of project sites. I am trying to make a new virtual field within the point layer that will reference the county that it is in. It needs to be able to stay in the attribute table of the point layer (and not creating a new joined layer) as it will be updated by adding more project sites/points in the future and am working in spatialite where this point shapefile has a join to a records table.

Does anyone know the expression I can use to do this?


texts: Connecting two coordinates with Polyline? Cartodb?

texts:
I have some trip data and I want to draw a line between the start and end coordinates.

I have two layers in CartoDB - one contains the start coordinates (with trip_id & cartodb_id). The other contains the corresponding end coordinates (with trip_id & cartodb_id).

As I'm new to CartoDB/SQL in general, how can I join these the two coordinates for the same trip?

I found this code somewhere but it hasn't been working for me thus far:

SELECT
  ST_MakeLine(
    the_geom_webmercator,
    (SELECT
      the_geom_webmercator
    FROM
      hubway_july2012_end
    WHERE trip_id = cartodb_id LIMIT 1
    )
  ) AS the_geom_webmercator
  FROM
  hubway_july2012_start
    AS c

(I've edited it a little bit, but I'm not sure if what I have done is right).


texts: Plotting GPS plot without uploading data

texts:
I'm not sure if this is possible or not, but I'm wondering if its possible to plot a GPS dataset without uploading the data first to CartoDB.

I've been analysing some data by uploading the postcode data to CartoDB and then using Javascript to control the mapping via the CSS and SQLApi which is really good. I can't however, figure out how to do this with several GPS plots, as I'd like to show a breadcrumb trail from data I have locally.


texts: Standard for determining accuracy of polygon interpolation technique?

texts: I'm looking into polygon interpolation, and was curious if there are any standards for determining the accuracy of polygon interpolation techniques. I've been looking through papers and lecturers slides and have found a reference (page 6) to removing sample points from the interpolation as a measure of accuracy. I've also found a paper (page 6) where they are using min, max, mean, and standard deviation as measures of accuracy. So, is there a standard out there for measure the accuracy of interpolation. Any published papers would be a big help.


texts: give importance weighting to sample points in an IDW-interpolation

texts:
I am trying to make an IDW-interpolation in arcgis 10.1.
In my sample file i have points with the same x,y and z values. According the arcgis-resource-site, the points are considered duplicates and have no effect on the output.

However, i want these points to have an effect:
the more often a certain value has been observed in a location, the more important that z-value should be.

Is there a way to circumvent this ignoring of duplicate sample points, or to use the number of duplicate points as a 'weight-value' to indicate relative importance of sample points?


texts: How to display previous epochs in of NLCD on GEE

texts:
GEE has the full National Land Cover Database and this has enabled a number of projects for me but I'm stumbling on a simple problem.

When I display the NLCD I can see that each point contains land cover values of all four epochs (1992, 2001, 2006 and 2011) but I'm not sure how to display them. As of now I basically just have the default display code but I'd like to show both 2011 and 2006.

var dataset = ee.ImageCollection('USGS/NLCD');
var landcover = dataset.select('landcover');
var landcoverVis = {
  min: 11.0,
  max: 95.0,
  palette: [
    '5475A8', 'ffffff', 'E8D1D1', 'E29E8C', 'ff0000', 'B50000', 'D2CDC0',
    '85C77E', '38814E', 'D4E7B0', 'AF963C', 'DCCA8F', 'FDE9AA', 'D1D182',
    'A3CC51', '82BA9E', 'FBF65D', 'CA9146', 'C8E6F8', '64B3D5'
  ],
};
Map.setCenter(-76.356, 35.286, 8);
Map.addLayer(landcover, landcoverVis, 'Landcover');


texts: Using multi iterate of ArcMap ModelBuilder

texts:
I want to use more than 2 iterates.
The data I have is point.gdb, polygon data.
There are about 2000 point.gdb.
There are about 1000 polygon data.
I want to save the points that exist in the polygon area for each polygon.
So I tried like this:

In ModelBuilder 1, I select the point data.gdb I have one by one (save as model1)

In ModelBuilder 2, iterates-feature selection is used to select polygon data one by one.

Using clip's function
input feature: read 1 (model1)
clip feature: read 2 times

Each item is set as a parameter

When I run it, it turns out that it cannot find the recursive function.
Where did I go wrong?


texts: Choosing target CRS value when doing geographic calculation?

texts:
In postgis, I use geography type, it support many geographic calculation, like st_dwithin(geography1, geography2, distance). What I need to do is just giving some points coordinates (lng, lat pair), on need to care about coord system.
Now I want to use java to achieve those functions.
I found this useful answer geotools-linestring-distance, @iant point out CoordinateReferenceSystem auto = CRS.decode("AUTO:42001,13.45,52.3");, after convert crs, I can get a correct distance.
But what "AUTO:42001" means? how can I choose a better target system?
Is there any documents or blogs to explain?


texts: pgRouting results to LINESTRING geometry type

texts:
I am doing some routing queries, using pgRouting's second version of pgr_trsp function.
What I can't do effectively yet is storing the results. That is, I want to able to store in a table:

either the actual result from pgRouting (not very useful anyway),or
geometry(LINESTRING) data type, derived from the
pgr_costresult (can be used in spatial queries).

Here is the bit of code I tried (and failed :P):

--Test function to convert pgr_costresult to geometry.
--Input arguments are the arguments of the pgr_trsp function, i.e. edge id's and relational position.
--TODO: do the same function work for arbitrary points.
CREATE OR REPLACE FUNCTION fm_pgr2geom(edge1 integer, pos1 float8, edge2 integer, pos2 float8)
    RETURNS geometry AS $$
 --We have to do a routing query. And declare a cursor for it
DECLARE resc CURSOR FOR
    SELECT * FROM pgr_trsp (
        'SELECT * FROM th_2po_4pgr',
        $1, $2, $3, $4, false, true);
    doline geometry[];
    temp_point geometry;
    geom geometry;
    temp_rec RECORD;
BEGIN
--Find startpoint and append to doline
    doline := array_append(
        doline,ST_Line_Interpolate_Point((SELECT map.geom_way FROM th_2po_4pgr map WHERE map.id = $1),$2));
--Append the rest of nodes
    FOR temp_rec IN SELECT * FROM pgr_trsp (
        'SELECT * FROM th_2po_4pgr',
        $1, $2, $3, $4, false, true) LOOP
            doline := array_append(
            doline, (SELECT nodes.geom FROM nodes WHERE nodes.id = temp_rec.id1));
    END LOOP;
--Append the endpoint
    doline := array_append(
        doline,ST_Line_Interpolate_Point((SELECT map.geom_way FROM th_2po_4pgr map WHERE map.id = $3),$4));
geom := ST_MakeLine(doline);
RETURN geom;
END;
$$ LANGUAGE PLPGSQL;

This piece of code is supposed to return a geometry (linestring) with all the points of the result. However the actual result is only a 2-point linestring.
Do you know any convenient way I can achieve either or both? I think second option is preferable, since I doubt postgresql tables can store complete pgr_costresults and Geometry types can be used for spatial analysis.


texts:
I want to update the 'Start_End' field in the following shapefile attribute table:

The attribute table consists of 6,134 records. For every unique value in the 'Id' field, I'd like to update each 'Start_End' record with a string, either 'Start' or 'End', depending on the value of 'OID_copy'. For example, for every unique value of 'Id' with the smallest value in the 'OID_copy' field, the 'Start_End' field should be populated with the word 'Start'. Conversely, for every unique value of 'Id' with the largest value in the 'OID_copy' field, the 'Start_End' field should be populated with the word 'End'.

Note: both the 'Id' and the 'OID_copy' fields are integer types, and the 'Start_End' field is string.

My initial thought was that a dictionary and an update cursor would do the trick, but so far I have been unsuccessful at updating the attribute table. My code doesn't throw an error, furthermore, the dictionary seems to have been created correctly. It looks like {1: [0, 1], 2: [2, 3], 3: [4, 5], ...}. However, the 'Start_End' field doesn't get updated.

# Start_End_UpdateAttribute.py
# purpose: For each unique 'Id' value, populate the 'Start_End' field in the valley centerline points shapefile attribute table
# with "start" and "end" using the value of the 'OID_copy' field as criteria

import arcpy
import traceback

arcpy.env.workspace = r'D:/GIS_data/Regional/Washington/ValleyBottom_UC_SF/ValleyBottom_UC_SF_centerlines/Wenatchee'
fc = r'D:/GIS_data/Regional/Washington/ValleyBottom_UC_SF/ValleyBottom_UC_SF_centerlines/Wenatchee' \
     r'/VC_ChiwawaRiver_split_verts_Copy.shp'
fields = ['Id', 'OID_copy', 'Start_End']
idDict = {}
cursor = arcpy.da.SearchCursor(fc, fields)

try:
    for row in cursor:
        id = row[0]
        oid = row[1]
        stend = row[2]
        idDict.setdefault(id, []).append(oid)
except:
    print 'An error occurred'
    traceback.print_exc()
del cursor
print idDict  # output is {key: [value, value], key: [value, value], etc}

cursor2 = arcpy.da.UpdateCursor(fc, fields)

try:
    for row in cursor2:
        for k, v in idDict.items():
            if v == min(v):
                row[2] = "Start"
                cursor2.updateRow(row)
            else:
                row[2] = "End"
    print 'Updated \'Start_End\' field'
except:
    print 'An error occurred - part 2'
    traceback.print_exc()
del cursor2

How should the second block of my code be revised to produce an attribute table that looks like this?:


texts: Convert JPG image to RGB Raster to Vector Data

texts: I am completing an analysis project for my ArcGIS class. I am analyzing the correlation between indigenous populations and food security in Guatemala. Since there is not a lot of information out there on this topic, I had to georeference .jpg maps I found online. I now have two RGB raster datasets: one shows the distribution of indigenous populations, and the other shows regional percentages of food scarcity. I would like to display the indigenous population raster layer as Quantities laid over the regional areas of food scarcity, therefore I have to convert it to vector data. I am having trouble finding a tool that can do this.


texts: Google Earth Engine VIIRS average radiance daily image

texts:
I'm using Google Earth Engine VIIRS data collection : "NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG",monthly average radiance data for VIIRS as:
var dataset = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG')
                  .filter(ee.Filter.date('2017-07-01', '2017-07-31')).map(clipToCol);
var nighttime = dataset.select('avg_rad');
var nighttimeVis = {min: 0.0, max: 60.0};

is it possible to use VIIRS daily data in Google Earth Engine?
If so how could be the code above modified to get the daily data or alternatively how can I access VIIRS radiance daily data?
Thanks a lot,
Gabriel


texts:
Can anyone explain how it is that the shp2pgsql loader for loading shapefiles into PostGIS goes about defining the schema for a table when importing a shapefile?

To further explain, shapefiles are comprised of geometry plus attributes of different data types and lengths: text, float, double, etc.  When running shp2pgsql with the -c option to create a table, it defines the data types on its own somehow, but it's not clear how this happens.  The only other option that looks related is -i, to use int4 for all integer fields.

Does it read into the file to find out the lengths of data values and consider those first?  Does it just assume basic things like text, integer, and float with default lengths?


texts: Where are the 'preloaded' images for the composer

texts:
I am trying to add a North arrow to a map in the composer in QGIS 2.0 on a Mac.

Googling tells me to simply add one from the preloaded images.  When I click the add image icon I get a dialogue box asking for a directory to find the images.  All the tutorials seem to assume that clicking on the add image button will bring up a bunch of handy images.  I am fairly new to QGIS and have not got this far with 1.8.

I can't find any docs specific to 2.0 and the print composer manual has not been updated yet.


texts:
Has anybody successfully used Zonal statistics in a loop in QGIS 3.8?

I have to go through several rasters using a vector file.

Batch processing seems to fail in QGIS 3.8 too - see Zonal statistics using Python console in QGIS gives AttributeError?


texts: Can custom map tiles be properly rendered via ArcObjects within my ArcGIS 10.1 Desktop AddIn?

texts:
Within ArcObjects for ArcGIS 10.1, we are attempting to render pre-generated map tiles as a
basemap (group) layer.  The tiles themselves are of the typical web mercator aux. spherical tiling scheme.

Under ArcMap 10.1 the tiles do not render properly:

and the ArcMap 10.1 UI is displaying the following drawing error: "Missing raster band wavelength properties"

[Expected] rendering result under ArcMap 10.0:

While this is related to this (my) question here, that was a holistic approach question.  Here I am specifically asking how to properly render the tiles in ArcMap 10.1.

We successfully accomplished the layer adding/rendering within 10.0 via the following code:

_tileLayer = new MyCustomTileLayer();    // implements ILayer among others..
_coverageBasemapLayer = new BasemapLayerClass();
_coverageGroupLayer = _coverageBasemapLayer as IGroupLayer;

// add the layer to the map
IMap map = ArcMap.Document.FocusMap;
map.AddLayer((ILayer)_coverageBasemapLayer);

...

// wire the AfterDraw Event
IActiveViewEvents_Event avEvent = ArcMap.Document.FocusMap as IActiveViewEvents_Event;
avEvent.AfterDraw += new IActiveViewEvents_AfterDrawEventHandler(avEvent_AfterDraw);

private void avEvent_AfterDraw(IDisplay Display, esriViewDrawPhase phase) {
    // force re-draw
    _tileLayer.Draw(esriDrawPhase.esriDPGeography, Display, null);
}

...

// effectively, within _tileLayer.Draw() here...
IRasterLayer rl = new RasterLayerClass();
// 'file' is the full path to the map tile PNG on the local file system, downloaded
// and/or cached
rl.CreateFromFilePath(file);

_coverageGroupLayer.Add(rl);

Presumably many things have changed within the renderer for 10.1.  However can anyone provide insight into the ArcObjects changes required?  Specifically how can we [properly] accomplish this within ArcObjects under 10.1?

Thank you in advance.


texts: Creating QGIS Installer for Windows with GUI customization

texts:
I'm trying to build a custom QGIS installer for Windows with a stripped down GUI. I customized the GUI using the customization function of QGis and saved it as a .ini file.

How can I include this custom GUI in my installer?

I found this very useful tutorial:

http://linfiniti.com/2011/05/building-custom-qgis-installers-for-windows/

and in the comments it also says that it's possible to include it using the customization tool. However, it doesn't exactly say how and where to save the customization file in the installer package.

Another possibility automatically load a custom GUI is explained here:

http://linfiniti.com/2011/07/batch-customizing-qgis/ Thought about exporting the registry    key customization.reg and to include it somewhere in the installer, but where?

So that would mean I would have to export the Windows registry key customization.reg, save it somewhere in the installation package and write a batch file that imports this registry key after the installation is finish.

Do you think that would be possible?

So that's what I have found out so far.

Does anybody have an idea what the easiest way would be to do this?


texts: Receiving higher values than should be possible in stackApply SUM for Frequency calculation in R

texts:
I am currently working on creating a frequency raster for a country based on multiple Rasters of the same day through following years and where the previous values are reclassified to 1's and 0's. The RasterStack is then SUMMED together to create the final frequency output. I am getting an output, but the problem is that the max value is too high.

For example, if I have 16 rasters in the stack, then the highest possible value should be 16 but instead I am getting a max value of 39. Below is the script I have created.

library(raster)

cloudstack <- stack(flsp) #Creation of RasterStack, 16 Rasters

m <- c(-Inf, 2, 0,  3, 3, 1,  4, Inf, 0) #Cloud = 1, Everything else = 0
reclassmatrix <- matrix(m, ncol=3, byrow=TRUE)
reclass <- reclassify(cloudstack, reclassmatrix) #Reclassify RasterStack

#SUM RasterStack
cloud_frequency <-stackApply(reclass,
             indices=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1), fun=sum, na.rm=TRUE)

crop_cloud <- crop(cloud_frequency, country)


texts:
when I try to create a workspace/layer/store following http://docs.geoserver.org/stable/en/user/restconfig/rest-config-examples/rest-config-examples-curl.html I am getting "System cannot find the file specified".... My curl is working fine with http.

What might be the solution?


texts: Error generating chart: No features contain non-null values of "system:time_start"

texts:
I know there are many posts abbout this but for some reason those solutions doesn't work for my code.

I'm trying to create timeseries for an imagecollection which was created using image expressiones.

I have tried to use

.copyProperties(im,['system:time_start']);

but didn't work and still got error message.

This is how I tried to keep the time properties in my code:

var FrostDB=clippedSAR.map(function(im){
   im = Frost(ee.Image(im),frostDamp,ksize);
   return im;
});

print(FrostDB);

var model=FrostDB.map(function(im){

  var nrpb = im.expression(
    'VH - VV / (VH + VV)', {
      'VH': im.select('VH_sum'),
      'VV': im.select('VV_sum')
    }).rename('NRPB').copyProperties(im,['system:time_start']);

  // adding NPRB Band to the original Image
  var withNRPB = im.addBands(nrpb);

  var result = withNRPB.expression(
      '2.572-(-0.05047*VH)+(0.17*VV)+(3.422*NRPB)', {
      'VH': withNRPB.select('VH_sum'),
      'VV': withNRPB.select('VV_sum'),
      'NRPB': withNRPB.select('NRPB').copyProperties(im,['system:time_start'])
    });

  // Only returning result band (other bands are not returned
  // you can add the result just like the nprb band was added if you want to return everything)
  return(result.copyProperties(im,['system:time_start']));
});

print('model:',model);

print(ui.Chart.image.series(model, geometry, ee.Reducer.mean(), 10));

My end goal:
To create time series without getting this error message.

edit: link to my code:
https://code.earthengine.google.com/4acde71b4b34093ac61cb76efe6a7339


texts: Why does Python only read the no data value from GeoTIFF?

texts:
I had a number of TIFF files with a country's map. They had different rows, columns or different projections etc. I needed to make them with same projection, spatial resolution and same number of rows and columns. So I used QGIS Raster->Projections-> Warp(Reproject) method and a vector boundary shapefile to get TIFF files with the intended parameters.

After doing this each file contained two bands. One band with the original data values and another band with values from 0 to 255. (I think it is for visualization, is it?)

In some of these files I had specified no data values and in others, I didn't. But when I'm reading these files in numpy arrays in Python using gdal, I'm only getting an array of intended size (53,54) but with  no data values! I used the following code:

from osgeo import gdal
from gdalconst import *
from scipy import stats
import numpy as np
import os, glob
import statsmodels.api as sm
import matplotlib.pyplot as plt

os.getcwd()
os.chdir(r"H:\...\Processed")
#load tiff files in numpy arrays
ds1 = gdal.Open("cases2001.tiff")
can01 = np.array(ds1.GetRasterBand(1).ReadAsArray()) #data type is float64
print can01.shape
ds2 = gdal.Open("cases2002.tiff")
can02 = np.array(ds2.GetRasterBand(1).ReadAsArray()) #data type is float 64
print can02.shape
#Load the RASTER CLIPPED REPROJECTED IMAGE AM I DOING SOME MISTAKE HERE?
ds3 = gdal.Open(r"H:\...\Annual 53-54\AOD2001.tiff")
AOD01 = np.array(ds3.GetRasterBand(1).ReadAsArray())
AOD01 =AOD01.astype('float64')  #Convertinf to float data type
print AOD01.shape
print AOD01

Why am I not getting the required values? And instead ending up with the "no data values"? For the files cases2001 and cases2002, which are the original files according to which I re-projected my other files are having the perfect, intended values. But for the third file here i.e. AOD2001 I'm getting an array of zeroes (I didn't change the default no data value). I've tried for both bands. But in QGIS "Identify Features" I'm getting the correct values for the 1st band. But not while reading ion Python.  I've already tried masking and a few other suggestions for similar problems (such asGetStatistics not refreshing etc) but I haven't found exact problem anywhere yet.

On a side note I need to plot these 2-D arrays against each other and perform regression. Which I'm obviously getting incorrect same values for each row as the AOD2001 file contains only zeroes. Another question, Can we compare (In my case perform regression on) two tiff files with same dimension, same spatial resolution and same end coordinates but different origins?

It's the first time I'm asking a question here, so I may have made mistakes in formatting etc., please point it out, I'd edit and take care in future. Also, please tell me if I should provide some other information or data.


texts: init.js:36 GET https://js.arcgis.com/3.28/esri/WebMap.js net::ERR_ABORTED 404

texts: I got above error in my browser developer tools window. This happened when I added esri/WebMap to my widget.js file of custom widget I was building.


texts: WFS bbox behavior for longitudes >= 360 degrees

texts:
I have a question regarding the OGC WFS specification (2.x and unfinished 3.x).

Let's say I have a WFS feature with longitude 0°.

If I make WFS query with bbox with minx = 359° and maxx = 361°, should it - according to the specification - return my feature?
In other words, should I expect, that 0° = 360°?

I'm also curious about the opposite functionality.
Given there is a service containing feature with longitude 360°, should this feature be returned through WFS query with bbox between -1° and +1°?

I'm not asking about particular technologies, but more about the expected behavior, according to the specification.


texts: Telling gdal_polygonize.py how to interpret colors / what does the `DN` field mean

texts:
I have lots of raster files I am converting to polys with gdal_polygonize.py. All these images have a consistent palette which I'd like to map to an attribute in the resulting polygons.

On most images, the output DN field is mapping consistently to the color values in the input rasters, however on one image it seems like the DN value is flipped.

I think it has something to do in which order the colors appear first in the image, because on that image, light orange appears first, then yellow and then orange (usually the order is orange, light orange, yellow and sometimes green).

What I've tried to fix it

The original input images are colored with a gradient from orange to green. I thought I could fix this by mapping these colors to a single band (I chose red) assigning the values #010000 through #040000 using convert ... -fill #010000 -solid -opaque #ff9900 ... and so on. This doesn't change the result though.
(in addition to 1) I'm already using gdal_edit.py to define the images bounding box, so additionally I've passed the flag -colorinterp_1 red to it and -band 1 to gdal_polygonize.py. This hasn't worked either.

Question

How do I tell gdal_polygonize.py how to interpret the colors in my raster files?

Images

This is the raster that fails to receive the correct values for the colors.

Orange = 3, light orange = 1, yellow = 2

This is an example raster that works as I expect.

Orange = 1, light orange = 2, yellow = 3

How the failing raster looks after mapping the palette to the red channel with convert (stretched with curves so it's not all black)


texts: Give non-overlapping parts of 2 rasters a value 0 (QGIS)

texts:
Is there a way in QGIS to give non-overlapping parts of 2 rasters a value of 0 in the new raster file?

I would like the data gaps in the velocity file (so where Ice mask is displayed on the chart) to be zero in a new raster file, so a new raster with the original value of velocity for pixels that have ice mask = 1 and velocity > 0, else 0.


texts: Save a temporary layer to a shapefile using PyQGIS

texts:
I'm attempting to save a polygon created in a scratch layer to a shapefile. The Python code below is a simplified example. Note the questions here and here provided some of the help, which I'm including in this post for future reference,
# Attach libraries
from qgis.core import * # attach main QGIS library
from qgis.utils import * # attach main python library
import os # attach operating system library

# Set the working directory
wd = "C:/test" # Set work directory
os.chdir(wd) # Change the directory

# Set a variable for the current project instance
  Prj = QgsProject().instance() # Object for current project

# Save the project to this file name
pnm = "Test.qgs" # Project file name
pnm = wd + "/" + pnm # Concat. with path
Prj.write(pnm) # Save the project

# Create an array [] object with the polygon vertices
vrtcs = []
vrtcs.append(QgsPointXY(396100,8969000))
vrtcs.append(QgsPointXY(396100,8973900))
vrtcs.append(QgsPointXY(397900,8973900))
vrtcs.append(QgsPointXY(397900,8969000))

# Create a polygon from the vertices
ply_01 = QgsGeometry.fromPolygonXY([vrtcs])

# Create a feature object then append the polygon into
ftr = QgsFeature()
ftr.setGeometry(ply_01)
print(ftr.geometry())

# Create a layer for the feature and add to the project
lyr = QgsVectorLayer('Polygon?crs=epsg:29194','Test',"memory")
Prj.addMapLayers([lyr])

# Make the layer editable, add the feature and save
lyr.startEditing()
lyr.addFeature(ftr)
lyr.commitChanges()

# Save as a shapefile
Fl_ou = 'Test.shp'
Fl_ou = wd + '/' + Fl_ou
writer = QgsVectorFileWriter(lry, Fl_ou,'urf-8','ESRI Shapefile')

When I execute the last line I get the following error  "NameError: name 'lry' is not defined".
However, when I enter type(lyr) into the console it recognizes the variable lyr as a <class 'qgis._core.QgsVectorLayer'> so I am unsure why I'm getting the not defined error.


texts: Many-to-Many relationship class table- automatically fill in?

texts:
Is there a way that the table of a many-to-many relationship class can automatically be filled in without me selecting the elements connected?
So far I am using this method: http://support.esri.com/en/knowledgebase/techarticles/detail/21708

But this takes lot of time, and I have huge amount of data, it's just not efficient to select every pole and connect it to every cable:P

PS: the cables and the poles are connected in a geometric network, can't arcgis figure out from that which one is connected to which one?


texts: How I can remove the gap between the mosaic image

texts:
I tried to mosaic the 2 Landsat images through
 (gdal_merge.py stacklayer.tif stacklayer_1.tif -o mosaaic.tiff)

following output come

How I can remove the gap although there projection is same

Before running the code both images looks like as above


texts: Showing startup-messages on QGIS splashscreen

texts:
During the startup of QGIS, there are status messages shown in the lower part of the splashscreen like "restoring loaded plugins".

I am using a startup.py feature from which I would like to inform the user which part of my startup script is executed at the moment.

Is it possible to show this information on the splashscreen?

Edit1:

As a workaround I managed to use my own splashscreen during startup:

from qgis.gui import *
from qgis.utils import *
from qgis.core import *
from PyQt4.QtGui import *
from qgis.PyQt.QtCore import QSettings, Qt
import time

template=QgsApplication.qgisSettingsDirPath() + "python/"
app=QgsApplication.instance()
splash_pix = QPixmap(template+'splashscreen.png')

splash = QSplashScreen(splash_pix, Qt.WindowStaysOnTopHint)
splash.setWindowFlags(Qt.WindowStaysOnTopHint | Qt.FramelessWindowHint)
splash.setEnabled(False)

splash.setMask(splash_pix.mask())

progressBar = QProgressBar(splash)
progressBar.setMaximum(10)
progressBar.setGeometry(0, splash_pix.height() - 20, splash_pix.width(), 10)

splash.show()

if QgsApplication.instance().findChild(QSplashScreen):
    QgsMessageLog.logMessage("ja", "gridseen", level=QgsMessageLog.INFO)
else:
    QgsMessageLog.logMessage("nein", "gridseen", level=QgsMessageLog.INFO)

splash.showMessage("
Grid Integration-Check!
", Qt.AlignBottom | Qt.AlignCenter, Qt.black)

for i in range(1, 11):
    progressBar.setValue(i)
    t = time.time()
    while time.time() < t + 0.1:
        app.processEvents()

time.sleep(2)
splash.close()

Therefore I put the splashscreen in my qgis-python folder (for example
https://github.com/webgeodatavore/qgis-splash-screens-birthday/raw/master/resized/qgis_version_2.18.png )

But this solutions is a bit of a quick and dirty workaround.

Isn't it possible to get access to the splashscreen that's created during the startup of the QGIS app? I tried to get access by using QgsApplication.instance().findChild(QSplashScreen) but couldn't get access to it.

https://github.com/qgis/QGIS/blob/7bd0285dfdef9456a5929a7b7031270ea0ee2601/src/app/main.cpp#L1286


texts: Adjusting space between symbols in legend using PyGIS and Print Layout

texts:
I have been trying to format a legend in order that it appears as a color bar under the map (in the print layout). I have been able to modify almost everything but I can't find a command to reduce the space between symbols to zero (mm).
Does anyone know how to do this?
The code until now is the following:
legend = QgsLayoutItemLegend(layout)
legend.model().setRootGroup(root)
layout.addLayoutItem(legend)
legend.attemptMove(QgsLayoutPoint(30, 200, QgsUnitTypes.LayoutMillimeters))
legend.setAutoUpdateModel(True)
legend.setItemRotation(270)
legend.setBackgroundEnabled(False)
legend.setFrameEnabled(False)
legend.setDrawRasterStroke(True)
legend.setSymbolHeight(11.0)
legend.setSymbolWidth(4.0)
legend.setEqualColumnWidth(0.0)
legend.setLineSpacing(0)
legend.setColumnSpace(0.0)
legend.setBoxSpace(0.0)


texts: Create Configurable Extension in ArcGIS Desktop

texts:
I am creating a configurable extension for ArcGIS desktop 10.1

I am not using Add-on , so please don't suggest add-on.

How can I disable some toolbars and commands, when the extension is disabled?
I want to disable users from adding toolbars from the "Customize" interface in ArcMap, when the extension is disabled?

Is there a way even to hide the toolbar from the "Customize" interface of GIS desktop when the extension is disabled?


texts: How to design a road with a certain slope

texts: Is there a way that you want to connect two different points of altitude with a constant slope of 7%  or another, to draw a line in the 3D model?


texts: Aligning two maps with different projections in QGIS

texts: I am trying to add historic borders to a modern blank world map. The files I have are all bitmap format and the historic borders are in separate regional maps (Africa, India, North America, Western Europe, etc.). I am quite new to this and am struggling considerably with accomplishing this. I was hoping to overlay the historic border maps as a transparent layer which would perfectly match the coastlines and rivers of the base map. I have attached two maps as a demonstration of what I mean (not in bitmap format as they were too large to do so).


texts: (openlayers) Adding styled polygons from predefined toolbar?

texts:
I'm looking for a solution where I could easily add differently styled polygons to a map?

Styles should be as "options" so that user can select e.g. blue or yellow polygon.

What different options there is to implement something like this? Is it to have two polygon buttons in toolbar? Or possibly some kind of attribute/style editing after polygons have been created?

The first one would be the most intuitive for end user and perhaps also easier for developer?

And what are the correct names for these ? From docs I can see that "controls" include "panel" which include "EditingToolbar controls" (http://dev.openlayers.org/docs/files/OpenLayers/Control/EditingToolbar-js.html).

How can I add another polygon button to toolbar that would use different custom styleset (or is stylemap in openlayers)?


texts: Loading Oracle Spatial Layers in standalone python qgis

texts:
I used below script to connect to postgres and load layers successfully.I want to aceive this using Oracle Spatial for both vector and raster loading.How can I achieve that?In other words what is the equivalent of postgis_utils for Oracle Spatial.From QGIS Desktop,I am able to load both raster and vector.This indicates QGIS has certain libraries to connect to Oracle and import GIS data.How can I use those libraries in standalone python script?

uri = QgsDataSourceURI()
uri.setConnection("localhost", "5432", "mydb", "postgres", "passcode")

db = postgis_utils.GeoDB(host="localhost", port=5432, dbname="v",user="postgres", passwd="passcode")
print db
tables =db.list_geotables()
print tables

render = QgsMapRenderer()

for t in tables :
   uri.setDataSource(str(t[1]),  str(t[0]), str(t[6]))
   uri.uri()
   vlayer = QgsVectorLayer(uri.uri(), str(t[0]), "postgres")
   QgsMapLayerRegistry.instance().addMapLayer(vlayer)
   lst = [ vlayer.id() ]
   render.setLayerSet(lst)
   rect = QgsRectangle(render.fullExtent())
   rect.scale(1.1)
   render.setExtent(rect)


texts: Getting ERROR 000840 from arcpy.SelectLayerByAttribute_management?

texts:
I am trying to select layer by attribute and for that, I am using arcpy.SelectLayerByAttribute_management function . My command is as follows:

arcpy.SelectLayerByAttribute_management(union, 'NEW_SELECTION', '"gridcode" =1 AND "area" >=25000000*0.4')

Here, union is the shapefile, from which I need to select layers by attribute.  My SQL query is "gridcode" =1 AND "area" >=25000000*0.4

I am getting the following error:

  ExecuteError: Failed to execute. Parameters are not valid. The value
  cannot be a feature class ERROR 000840: The value is not a Raster
  Layer. ERROR 000840: The value is not a Mosaic Layer. Failed to
  execute (SelectLayerByAttribute).`


texts:
I am using the KML tools pluggin to import kmz files but I get:
"Failed to open file
Failed to open file
Execution failed after 0.01 seconds"
What could be causing this error?


texts: Merge overlapping polygons based on attribute

texts:
My objective is very similar to these questions: here and here. The only real difference is I want to merge overlapping polygons if they are of the same class (an attribute).
For example, if polygons of different classes A and B overlap, they should be left as is, while if both are A or both are B, they should merged/dissolved.
                                            geometry  class
0  POLYGON ((233572.138 2710742.125, 233576.808 2...    A
1  POLYGON ((233539.259 2710809.727, 233543.978 2...    A
2  POLYGON ((233509.510 2710826.481, 233514.245 2...    A
3  POLYGON ((233571.415 2710766.758, 233576.102 2...    B


texts: Geoserver multiple nodata values ECW in alpha band

texts: I want to combine multiple geotiffs to one ECW in FME and then upload the ECW in Geoserver. The responding nodata is stored in the alpha band. However with ECW compression (even when choosing 0 for compression) the 0 value suddenly becomes 1 or more so in Geoserver the background edges are only partially transparent (unlike in QGIS where the entire alpha band can be declared as nodata). Is there a way to declare all values in the alpha band as nodata in Geoserver (or a range) or is there another solution to this in FME or Geoserver?


texts: ModelBuilder Run Time Errors Using ArcGIS 10.0

texts:
I recently susccesfully ran a watershed delineation model using both ArcHydro Tools for Arc GIS 10.0 and Model Builder in ArcGIS 10.0 on Windows 7 system (Dell Inspiron N4110 Intel Core TM i3 CPU @ 2.30 GHz 64 Bit.)

However when I try to run the same model on an XP computer (Intel Pentium 4 CPU @ 2.53 GHz 32 Bit) I received several run time errors such as:

Error HRESULT E_FAIL has been returned from a call to a COM component. Fail to execute (Catchment Polygon Processing)" and also other like
Drainage Point.. All the inputs are not current
Adjoint Catchment All the inputs are not current

I try also to build the same model from the scratch in this later system and I received the same type of run time errors.  I'd be pleased to receive any clue to fix this problem.


texts: Manually set QGIS vertex coordinate in lat-long format

texts:
I'm new to QGIS (and GIS in general)

I'm creating a map and would like to manually enter lat-long coordinates in decimal format.

As you can see, the vertex editor uses a different coordinate system (I assume it's metres) but I'd like to set coordinates in decimal format, just like the bottom right of the photo

How can I do this?


texts: Mapcache: how to calculate the x,y position of a tile?

texts:
I use MapCache of Mapserver to generate tiles of a tif.

Given a tile (located on disk as z/y/x structure), how can I calculate it's geo-position?

Example:  8/50/100.png ->  ( 46.11, 11.05 ) (in wgs84)


texts: Summarizing table columns in a dynamic text expression

texts:
I have some data in a table that I would like to perform some statistics and have that placed in a text box on my map Layout. While I could do this by hand, I'd like to do it dynamically so when I update the data later on I don't have to edit my text boxes.

My data looks something like this:

| County | POP2019 | Deaths |
|--------|---------|--------|
| Benton | 1000    | 2      |
| King   | 100     | 1      |

I want to do a statewide death rate calculation like this

((Sum of Deaths / Sum of Population) * 100) Round to 5 decimals

Round(((Sum($feature.Deaths) / Sum($feature.POP2019)) * 100), 5)

But that dynamic text expression doesn't sum up the column! It prints out each and every row, and does not sum them up.


texts: How do I get road intersection symbology to display properly?

texts: so, I'm finding that at intersections of roads, I am unable to get correct display by making changes in the road symbols themselves in Arc 10. Sometimes one road will display on top of the intersecting road, and sometimes vice versa--I haven't found any apparent pattern to the display. There must be another way to do this. ???


texts: Computing mean direction of river network in R

texts:
I have a shapefile of a section of river network (polyline) which I have read into R using the sp package to have a SpatialLinesDataFrame. I also have access to a terrain raster for the region in which the river section lies. I wish to find the mean direction of the river network (analogous to ArcGIS's Linear Directional Mean) but I want to keep it in R.

I can't see the function I want in the sp or rgdal or raster packages. Is there a package with the appropriate function in it?

PS.
This (Is there an R alternative for ArcMap's linear directional mean tool?) question shows the brute force method using a table with x_start, x_end, y_start, y_end. Is there anything better than this? If not, how do I extract such a table from a SpatialLinesDataFrame?


texts: Consequences to QGIS of changing PostGIS table to 3D?

texts:
In QGIS add vector layer creates a layer with 3D geometry.  All my postGIS tables have been created with 2D so I can not save them to PostGIS.

Are there any problems with converting my existing 2D tables to 3D to get around this?

I am aware that QGIS 3 will allow one to drop the z dimension but this question is explicitly about converting tables to 3D.


texts:
I have a QGIS project which contains a composer. Of course if I go into the GUI I can export it as a pdf etc. For my make script though I want to have this done via a python script. From reading online here is my attempt:

#!/usr/bin/env python

import sys
from qgis.gui import QgsMapCanvas, QgsLayerTreeMapCanvasBridge
from PyQt4.QtCore import QFileInfo
from PyQt4.QtXml import QDomDocument
from qgis.core import *

qgs = QgsApplication(sys.argv, True)
QgsApplication.setPrefixPath("/usr", True)

qgs.initQgis()

def make_pdf():
  canvas = QgsMapCanvas()
  bridge = QgsLayerTreeMapCanvasBridge(
      QgsProject.instance().layerTreeRoot(), canvas)
  bridge.setCanvasLayers()
  QgsProject.instance().read(QFileInfo('../board.qgs'))

  composition = QgsComposition(canvas.mapSettings())
  map_item = composition.getComposerItemById('board36x48')
  map_item.setMapCanvas(canvas)
  map_item.zoomToExtent(canvas.extent())
  composition.refreshItems()
  composition.exportAsPDF('generated/board.pdf')
  QgsProject.instance().clear()

make_pdf()

This fails, specifically:

i:./assets/generate_board
QGraphicsScene::addItem: item has already been added to this scene
Traceback (most recent call last):
  File "./assets/generate_board", line 30, in <module>
    make_pdf()
  File "./assets/generate_board", line 24, in make_pdf
    map_item.setMapCanvas(canvas)
AttributeError: 'NoneType' object has no attribute 'setMapCanvas'

I based this off of Save Print/Map QGIS composer view as PNG/PDF using Python (without changing anything in visible layout)? but modified it to simplify as much as possible. It seems like I broke something in the process... I don't have this template file. I simple have a board.qgs project file with a composer board36x48. Any ideas?


texts: Local variable 'value' referenced before assignment: qgis function editor

texts:
I customized a function to edit my attribute table. It was working well but now when I used it again my QGIS crash or I receive a bad answer.

The answer of QGIS:

My code:

from qgis.core import *

@qgsfunction(args="auto", group='Custom')
def change_name_field(field, feature, parent):
    if 'agropecuaria' in str(field):
        value1 = 'Agroforestry'
    elif 'agua' in str(field):
        value1 = 'Water'
    elif 'cabruca' in str(field):
        value1 = 'Agroforestry'
    elif 'capoeira' in str(field):
        value1 = 'Secondary Forest'
    elif 'eucalipto2limao' in str(field):
        value1 = 'Silviculture'
    elif 'eucalipto3limao' in str(field):
        value1 = 'Silviculture'
    elif 'eucalipto2mad' in str(field):
        value1 = 'Silviculture'
    elif 'eucalipto4madesc' in str(field):
        value1 = 'Silviculture'
    elif 'mangue' in str(field):
        value1 = 'Mangrove'
    elif 'mata' in str(field):
        value1 = 'Forest'
    elif 'nuvem' in str(field):
        value1 = 'Cloud'
    elif 'solo arenoso' in str(field):
        value1 = 'Sand Soil'
    elif 'solo exposto' in str(field):
        value1 = 'Bare Soil'
    elif 'solo umido' in str(field):
        value1 = 'Soil Moist'
    elif 'sombra' in str(field):
        value1 = 'Shadow'
    elif 'Agroforestry' in str(field):
        value1 = 'Shadow'
    return value1


texts: Dissolve spatio temporal data with GeoPandas without missing the temporal part

texts:
I am wondering if there's a clever way to use the dissolve() method of GeoPandas with a spatio temporal geodataframe without losing the temporal part.
I have a dataset where for different days of the year the same geometries could repeat. I need to aggregate (dissolve) these geometries by ID but I need to do it for each day.
My geodataframe looks like this:
geometry DATE ID VOL
POINT 2020-01-01 1 100
POINT 2020-01-01 1 50
POINT 2020-01-01 2 50
POINT 2020-01-01 2 50
POINT 2020-01-02 1 70
POINT 2020-01-02 1 40
POINT 2020-01-02 2 30
POINT 2020-01-02 2 70
...

My desired output should be:
geometry DATE ID VOL
MULTIPOINT 2020-01-01 1 150
MULTIPOINT 2020-01-01 2 100
MULTIPOINT 2020-01-02 1 110
MULTIPOINT 2020-01-02 2 100

With gdf.dissolve(by='ID', aggfunc='sum') the result is missing the temporal information, and I understand why, but, is there a way I can aggregate preserving it?
geometry ID VOL
MULTIPOINT 1 260
MULTIPOINT 2 200

Here is the code in order to replicate the test:
import geopandas as gpd
from shapely.geometry import Point

d = {
    'geometry': [Point(1,2), Point(2,1), Point(3,4), Point(4,3),
                 Point(1,2), Point(2,1), Point(3,4), Point(4,3)],
    'DATE': ['2020-01-01', '2020-01-01', '2020-01-01', '2020-01-01',
             '2020-01-02', '2020-01-02', '2020-01-02', '2020-01-02'],
    'ID': [1, 1, 2, 2, 1, 1, 2, 2],
    'VOL': [100, 50, 50, 50, 70, 40, 30, 70 ]
}

gdf = gpd.GeoDataFrame(d, crs="EPSG:4326")
dissolved = gdf.dissolve(by='ID', aggfunc='sum')


texts: QGIS Time Manager 2.2.1 does not display any points

texts:
I have a large data set collected by a GPS collar on a cow that I am trying to animate with the Time Manager plugin to show the cow's movement throughout a pasture. I have the data saved in a .shp file that displays correctly. When I turn on the TM plugin and go through the settings then click "OK" all my points disappear from the canvas. I click "play" and the plugin appears to progress through time, but no points display. I noticed when I try to open the attribute table while the plugin is on nothing displays, but when I switch it off everything reappears.

Attached are some screenshots of my settings and what I am seeing:
http://i.imgur.com/RVVGOGt.png Settings

http://i.imgur.com/HtHzDjl.png TM on

http://i.imgur.com/uYHesps.png TM off

I have double checked time format and it is in YYYY/MM/DD HH:MM format. I can see the time stamp in the corner of the map progressing through but no points display at all.


texts: Automatic folder selection for processing algorithms in graphical modeller

texts:
Background

I have 19 layers which ALWAYS have the same file name, but are located in varying folders. These layers are all inputs to several big models.

Now, instead of defining each algorithm input manually, I want my model to pick them up automatically by only providing a folder location.

Basically I create one algorithm for each layer which has the file name hard-coded, it takes a an input folder location from the model and outputs a layer (code for an example algo at the end). But even now, I have to tell each algorithm manually to use the model's input folder parameter, schematically (in reality there's 16 more algos needing the same folder input)

Question

Is there a way that I can tell the "Get OUT_Distribution.." algorithms to use the first folder picker parameter it finds if the environment is a model? I expect it to be available over QgsProcessingAlgorithm.initAlgorithm's configuration argument, but there's little guidance on how to use that.

Or maybe this is overly complex for what I want to achieve?! @Joseph had a nice solution via simple Python scripts, but I've NEVER seen that notation before, nor can I find any documentation on how one could add inputs using that method.

Example code

Example code for a custom algo returning a layer from hard-coded file name and input folder location (just to show the gist of it, not necessarily working):

class ExampleAlgo(QgsProcessingAlgorithm):

    IN_FOLDER = 'IN_FOLDER'
    OUT_LAYER = 'OUT_DistributionCableEntries'

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFile(
                self.IN_FOLDER,
                'Input folder',
                behavior=QgsProcessingParameterFile.Folder,
                fileFilter='Alle Dateien (*.*)',
                defaultValue=None
            )
        )

        self.addParameter(
            QgsProcessingParameterVectorDestination(
                self.OUT_LAYER,
                'OUT_DistributionCableEntries',
                type=QgsProcessing.TypeVectorAnyGeometry,
                createByDefault=True
            )
        )

    def processAlgorithm(self, parameters, context, model_feedback: QgsProcessingFeedback):

        input_folder = self.parameterAsString(
            parameters,
            self.IN_FOLDER,
            context
        )

        out_layer = self.parameterAsVectorLayer(
            parameters,
            self.OUT_LAYER,
            context
        )

        output_path = path.join(input_folder, self.OUT_LAYER) + '.shp'

        layer_out = QgsVectorLayer(output_path, self.OUT_LAYER, 'ogr')

        return {self.OUT_LAYER: layer_out}

    def name(self):
        return f'get_{self.OUT_LAYER}'

    def displayName(self):
        return f'Get {self.OUT_LAYER}'

    def createInstance(self):
        return D4SAlgoGetFileDistributionCableEntries()


texts: Qfieldsync : Using multiple field devices

texts:
I need to use Qfield to collect data from multiple field workers. Is this possible with Qfieldsync and if so what would the workflow be.

If I sync the first device with the desktop project what will be the result of synching the second device; will it overwrite the data from the first device?


texts:
When I run gdalinfo (v2.2.2) on a geotiff that I created from a pdf using gdal_translate, I don't get the pixel size as shown in the docs. Does the pixel size have to be explicitly stored in the tif in order for gdalinfo to report it? If so, can gdal_translate be coerced to store the value when the tif is written?

>gdalinfo NM_Canada_Ojitos_20170216.tif
Driver: GTiff/GeoTIFF
Files: NM_Canada_Ojitos_20170216.tif
Size is 3412, 4350
Coordinate System is:
PROJCS["NAD83 / UTM zone 13N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-105],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","26913"]]
GeoTransform =
  319569.9072199312, 4.06320686118055, -0.08151911420909382
  4042818.402645736, -0.08151911420909382, -4.06320686118055
Metadata:
  AREA_OR_POINT=Area
  AUTHOR=USGS National Geospatial Technical Operations Center
  CREATION_DATE=D:20170216103259Z
  CREATOR=ESRI ArcSOC 10.0.2.3200
  KEYWORDS=Topographic, Transportation, Hydrography, Orthoimage, U.S. National Grid, imageryBaseMapsEarthCover, Imagery and Base Ma
s, Geographic Names Information System
  NEATLINE=POLYGON ((332137.201278231 4041102.99001007,331856.525171518 4027113.079869,320528.939318619 4027340.34242206,320809.615
25332 4041330.25256312,332137.201278231 4041102.99001007))
  SUBJECT=This image map depicts geographic features on the surface of the earth.  It was created to provide a representation of ac
essible geospatial data which is readily available to enhance the capability of Federal, State, and local emergency responders for
omeland security efforts.  This image map is generated from selected National Map data holdings and other cartographic data.
  TITLE=USGS 7.5-minute image map for Canada Ojitos, New Mexico
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  319569.907, 4042818.403) (107d 0'53.88"W, 36d30'49.40"N)
Lower Left  (  319215.299, 4025143.453) (107d 0'53.29"W, 36d21'15.87"N)
Upper Right (  333433.569, 4042540.259) (106d51'36.58"W, 36d30'49.43"N)
Lower Right (  333078.961, 4024865.310) (106d51'37.13"W, 36d21'15.87"N)
Center      (  326324.434, 4033841.856) (106d56'15.22"W, 36d26' 2.73"N)
Band 1 Block=3412x1 Type=Byte, ColorInterp=Red
Band 2 Block=3412x1 Type=Byte, ColorInterp=Green
Band 3 Block=3412x1 Type=Byte, ColorInterp=Blue


texts: Creating/Deleting fields (with expression) for a shapefile in a batch script

texts:
I am new to GDAL/Scripting in general.
For a study project, I have to create two new fields to a shapefile/vector file using a certain expression. I also need to delete two fields.
Normally, I would just use the Field Calculator in QGIS, choose "Create a new field" and delete the two other old fields. But I have to automate the process, and as there is no command for the field calculator in GDAL; I do not know what to do. Surely, you could do something with ogr2ogr, but unfortunately I do not find any clue in the documentation there either.
How do I add the expression, which creates the new columns?
It is the following:
case
when  "COVERTYPE" =  'AGR' then 1
when  "COVERTYPE" =  'BAR' then 2
when  "COVERTYPE" =  'CON' then 3
when  "COVERTYPE" =  'HDW' then 4
when  "COVERTYPE" =   'HEB' then 5
when  "COVERTYPE" =  'MIX' then 6
when  "COVERTYPE" =  'SHB' then 7
when  "COVERTYPE" =  'URB' then 8
when  "COVERTYPE" =  'WAT' then 9
end

Also, the output field type of the two new columns have to be "Whole number (integer 64 bit)".
How do I formulate that for the script?
I'm using QGIS 3.16.3.
UPDATE: I tried it now with the following commands:
ogr2ogr -append -sql "alter table vegetation_two_columns ADD column NVCS_CLASS_NEW INTEGER(64)" "Wildfireanalysis_Santa-Cruz\Vegetation\Modified_columns\vegetation_modified.dbf" "Wildfireanalysis_Santa-Cruz\Vegetation\Irrelevant_data_deleted\vegetation_two_columns.dbf"
and
ogr2ogr -update -append -overwrite -sql "alter table vegetation_two_columns ADD column NVCS_CLASS_NEW INTEGER(64)"  "Wildfireanalysis_Santa-Cruz\Vegetation\Modified_columns\vegetation_modified.dbf" "Wildfireanalysis_Santa-Cruz\Vegetation\Irrelevant_data_deleted\vegetation_two_columns.dbf"
While executing both of the commands, thw following error message appears: "ERROR 6: CreateField : unsupported operation on a read-only datasource."


texts: Creating density map from point data in ArcGIS Desktop?

texts:
I have point data, where one point represents one 1 storm water outfall, looks like this:

(point data is across all of Long Island, New York).

I am having a bit of difficulty finding the best tool to create a raster density map. Essentially I want to use a density layer in a weighted overlay model (an area with high density of storm water outfalls would be less suitable...). I do not think I need to use Hot Spot Analysis because I am only interested in high v low density and not if they are actually statistically different from each other. I also have tried point density tool but I am having trouble figuring out the inputs (like radius length and cell size, etc).


texts: Convert Nasa's Blue Marble into a cubic projection

texts:
I'm looking for a free tool to convert the cylindrical projection map of the Earth's surface that you can get from Nasa (under the Blue Marble project):
http://visibleearth.nasa.gov/view_cat.php?categoryID=1484

to a cubic projection like so:
http://www.progonos.com/furuti/MapProj/Normal/ProjPoly/Foldout/Cube/cube.html

I'd appreciate any assistance...


texts: Supporting group editing of same shapefile with open source tools?

texts:
We are a small (~5) group working on a single, shared shapefile. The current workflow is that each one of us does the editing on a designated attribute field, on his local copy, then exports the attributes table as a CSV file, which will be joined with shapefiles on all the other copies.

Obviously it is quite inefficient and, when comes to adding a feature or multiple people working on the same attribute, very error-prone.

How can we manage this situation?


texts:
I am see differences between tables when I view them in spatialite-gui and in QGIS. I have created a table, added columns, and updated the columns in spatialite-gui. Within spatialite-gui I can view all the added columns, as well as the updated fields. However, when I loaded the same table into QGIS 2.2 and open the attribute table, I do not see the columns I have added. I am able to see all the columns in QspatiaLite, but not the attribute table in QGIS (see graphic below).

I am using a Windows 7 64-bit computer with QGIS 2.2 from OpenGeoSuites4W, and spatialite-gui 1.7.1

The following are the queries I have written to create, add columns, and update the table.

 Creating table with three columns, then adding geometry, and finally inserting geometry into the table from another using a union statement. This was first performed in spatialite- gui 1.7.1

 CREATE TABLE smc_bldg2d (
     gid INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT,
     city TEXT NULL, state TEXT NULL, country TEXT NULL);

 SELECT AddGeometryColumn('smc_bldg2d', 'geom', 4326, 'MULTIPOLYGON', 'XY');

 INSERT INTO smc_bldg2d
     (city, state, country, geom)
 SELECT f.city AS city, STATE AS state, COUNTRY AS country, st_unaryunion(st_collect(f.Geometry)) AS geom
 FROM smc_bldg2d_usa_poly AS f
 GROUP BY city
 ORDER BY city;

 ALTER TABLE smc_bldg2d ADD COLUMN type TEXT;

 UPDATE smc_bldg2d
 SET type = 'vector'
 WHERE gid > 0;


texts: Calculating max shardSize and fileDimensions to avoid tiling large raster exports in Earth Engine

texts:
I am exporting a very large raster in Google Earth Engine to my Google Drive.  I am trying to keep my exported raster to 1 file to avoid additional processing to merge the tiles (e.g. GDAL merge in QGIS).  I know you can specify larger shardSize and fileDimensions to increase the size of your exported raster tiles.

For example, based on this post, I've been setting my shardSize to 1024 and playing around with my fileDimensions by trial and error.  So I arbitrarily set my fileDimensions to 1024*17 = 17408 but return an error on the following...

  Export.image.toDrive({
    image: large_raster
    scale: 30,
    maxPixels: 1e13,
    region: aoi,
    shardSize: 1024, // must be no greater than 1024; larger shardSize --> larger tile
    fileDimensions: 17408 // must be a multiple of 1024 and no greater than 131072
  });

Error: File dimensions are too large for the number of bands and data type: the total file size requested is 32728154112 bytes before compression (108 bytes per pixel), but must be less than 17179869184 bytes.

...then I will adjust my fileDimensions to the next lowest multiple of 1024 (e.g. 1024*16=16384), try to export, and see if I return an error.

But, I want to calculate my max fileDimensions from the get-go without this trial and error.  Does anyone know how to do that?


texts: Issue converting projected coordinates to wgs1984 in Postgis

texts:
I'm trying to convert a local PCS (city of Buenos Aires) to wgs1984. Since the PCS is not included in Postgis, I added it to spatial_ref_sys using spatialreference.org (I compared the arguments with a .prj I know is right and they match).
My coordinates for the PCS are in WKT, so to make my geometry column in wgs I used:

update <my_table>
set geom = ST_Transform(ST_SetSRID(ST_GeomFromText(<wkt_column>), <PCS_number>), 4326);

Which works at first glance, but the points are off by approx. 50 meters in relation with my other wgs tables.

When I make the transformation using the Project tool on Arcgis, it works perfect. The difference I see is that Arcgis offers me to use different 'Geographic transformation methods' to make the transformation. The one I choose, and it works fine, is 'Campo_Inchauspe_To_Posgar_2007_1 + Posgar_2007_To_Wgs_1984_1'. I suspect that this method is the key to solve this, but in that case I don't have a clue on how to include that method in the Postgis transformation.

My .prj is this:

PROJCS["Argentina_GKBsAs",
GEOGCS["GCS_Campo_Inchauspe",
DATUM["D_Campo_Inchauspe",
SPHEROID["International_1924",6378388.0,297.0]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",100000.0],
PARAMETER["False_Northing",100000.0],
PARAMETER["Central_Meridian",-58.4627],
PARAMETER["Scale_Factor",0.999998],
PARAMETER["Latitude_Of_Origin",-34.6297166],
UNIT["Meter",1.0]]

The projection number in spatialreference.org is 8333, and here you can see it in wkt as html for comparison:

  PROJCS["Gauss Krugger BA",
      GEOGCS["GCS_Campo_Inchauspe",
          DATUM["D_Campo_Inchauspe",
              SPHEROID["International_1924",6378388.0,297.0]],
          PRIMEM["Greenwich",0.0],
          UNIT["Degree",0.0174532925199433]],
      PROJECTION["Transverse_Mercator"],
      PARAMETER["False_Easting",100000.0],
      PARAMETER["False_Northing",100000.0],
      PARAMETER["Central_Meridian",-58.4627],
      PARAMETER["Scale_Factor",1.0],
      PARAMETER["Latitude_Of_Origin",-34.6297166],
      UNIT["Meter",1.0]]

Here is the Postgis insert statement:

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 98333, 'sr-org', 8333, '+proj=tmerc +lat_0=-34.6297166 +lon_0=-58.4627 +k=1 +x_0=100000 +y_0=100000 +ellps=intl +units=m +no_defs ', 'PROJCS["Gauss Krugger BA",GEOGCS["GCS_Campo_Inchauspe",DATUM["D_Campo_Inchauspe",SPHEROID["International_1924",6378388.0,297.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",100000.0],PARAMETER["False_Northing",100000.0],PARAMETER["Central_Meridian",-58.4627],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",-34.6297166],UNIT["Meter",1.0]]');

So, I have one point in wgs:

-58.47532313, -34.63847446

And I make the geometry (Point, 4326) in one table

Then I take the same point in the PCS (8333):

POINT(98897.6171875 98978.84375)

And I make the geometry (Point, 4326) out of the following transformation:

update test_t
set geom = ST_Transform(ST_SetSRID(ST_GeomFromText(coord), 98333), 4326);

I get this when I pick the point in Qgis:

Being the green point the correct one (the one originally in wgs), and the red the misplaced one (the transformed).


texts: I have to build a model for my class midterm which requires us to take a .csv file and convert it to a table, then take the table and create an XY event layer which will be used to create a new feature class.  From that new feature class, I am to break the trees into separate feature classes by tree species which will be placed into a feature dataset.  When I run the model as built with the first .csv file for Lake Shore Trees, it runs as it should and I get exactly what I was expecting.  The problem I keep running into is when I go to run the model again and this time change the allowable parameters for the next .csv file to be run, I keep getting "Error 999999" and the only other additional information is "failed to execute Feature Class to Feature Class."  When I look in my catalog tab, I see that the model has added only one feature to the feature NorthPark dataset, when there should be a total of three features in the dataset.  Now if I build a completely new model and use NorthPark as my initial input .csv file for the model, it runs as expected and all features are added to the dataset, but I get the same error again when I go and change the input parameters for the LakeShore .csv file.  Anyone have an explanation as to why this might be happening?


texts: OziExplorer maps not rendered correctly in QGIS

texts:
When I open up some OziExplorer (.ozf2) maps in QGIS they are not rendered correctly in that they appear reversed or in-versed. Text on the map is not legible and looks to be in-versed but this is not the case (looking via a mirror confirms this). Find this very strange and looking for some help.

Maps from http://maps.festy.org/


texts: Is there a fix for "sqlite3_exe(COMMIT) failed: database is locked" QGIS 3.6 Python 3.7

texts:
I am getting an error "ogr error committing transaction: sqlite3_exe(COMMIT) failed: database is locked" when trying to edit a feature on a geopackage layer using pyqgis.    It is a QGIS plugin.   Everything works fine on my hard drive, the problem occurs when I load the file to a server.

I have found https://github.com/r-spatial/sf/issues/628   which seems similar.
This also - How to update attributes on spatialite table programmatically: SQLite Error

def save_point(self):
    try:

        node_lyr_name = 'points'

        node_point = QgsPoint(-83.433, 44.466, 725.000)
        point_id = "2001"

        node_lyr = QgsProject.instance().mapLayersByName(node_lyr_name)[0]

        features = node_lyr.getFeatures()
        node_lyr.startEditing()
        for feature1 in features:
            if feature1['ID'] == point_id:
                feature1.setGeometry(node_point)
                feature1['Name'] = self.dlg.le_yourname.text()
                node_lyr.updateFeature(feature1)
                break

        node_lyr.commitChanges()
        node_lyr.triggerRepaint()

    except Exception as inst:
            self.popMsg(str(inst) + "\n Exception #1561 for row in save_point(): \n" )
            print(type(inst))
            print(inst.args)
            print(inst, ' #1604')

Is there anything different I need to do when working with a geopackage?

[EDIT]  the below code also fails with the same error.  :(

vectorLayer = iface.activeLayer()

features = vectorLayer.getFeatures()
fid = int()
    for feature in features:
    if feature['ID'] == "2016":
        fid= feature.id()
        break

field = vectorLayer.fields().indexFromName('Name')
print(f'Field {field}')
with edit(vectorLayer):
    vectorLayer.changeAttributeValue(fid, field, 'Hello')


texts:
I have received a dataset projected to Lambert Conformal Conic though it appears to be a custom projection.  The layer properties:

  Projected Coordinate System:   Lambert_Conformal_Conic

  Projection:    Lambert_Conformal_Conic

  false_easting: 600000.00000000

  false_northing:    2200000.00000000

  central_meridian:  0.00000000

  standard_parallel_1:   45.89891896

  standard_parallel_2:   47.69601450

  latitude_of_origin:    46.80000000

  Linear Unit:   Meter

  Geographic Coordinate System: GCS_unnamed ellipse

  Datum: D_unknown

  Prime Meridian: Paris

  Angular Unit: Degree

This is causing me issues when trying to integrate other datasets in WGS84 and reprojecting.

Edit: here is the processing error in ArcGIS 10.4:

And here is the output if I proceed with the projection (data should overlay):

Is there a solution to this other than trying to find out the original source for the dataset?


texts: Extract median value for polygons from multiresolution raster data, using R (velox, raster -package)

texts:
I am trying to extract median values from rasters for 1043 different polygons (plots).

Raster resolution varies from 1 to 30 meters and the plot size is r = 6 meters.  As you can see, at coarser resolution data the polygon doesn't always contain any pixel center.

The polygons area in a SpatialPolygonsDataFrame and I am trying to use the veloxpackage for the extraction. For the smaller resolution it works beautifully;

v1$extract(spdf,fun = median) # v1 = veloxraster, 1 m resolution

returns the median value for the plots. However;

v30$extract(spdf,fun = median) # v30 = veloxraster, 30 m resolution

returns the median for those plots which are divided to multiple pixels (or pixel centers), and returns NA for those polygons, which do not coincide with any pixel center coordinates.

The raster package has an option small = TRUE, where the single pixel value is returned. According to the packages' documentation:

  Logical. If TRUE and y represents points and a buffer argument is
  used, the function always return a number, also when the buffer does
  not include the center of a single cell. The value of the cell in
  which the point falls is returned if no cell center is within the
  buffer. If y represents polygons, a value is also returned for
  relatively small polygons (e.g. those smaller than a single cell of
  the Raster* object), or polygons with an odd shape, for which
  otherwise no values are returned because they do not cover any raster
  cell centers.

What would be an easy way to achieve that with velox? The extraction speed difference is huge (1.5 s / raster vs. 360-380 s / raster), so I'd like to perform these operation with velox if by any means possible.

# Reproducible example:
library(velox)
library(raster)

## Make VeloxRaster
mat <- matrix(1:100, 10, 10)
extent <- c(0,1,0,1)
vx <- velox(mat, extent=extent, res=c(0.1,0.1), crs="+proj=longlat +datum=WGS84 +no_defs")

## Make SpatialPolygonsDataFrame
library(sp)
library(rgeos)
set.seed(0)
coords <- cbind(runif(10, extent[1], extent[2]), runif(10, extent[3], extent[4]))
sp <- SpatialPoints(coords)

# Default example
# https://cran.r-project.org/web/packages/velox/README.html
spol_norm <- gBuffer(sp, width=0.2, byid=TRUE)
spdf_norm <- SpatialPolygonsDataFrame(spol_norm, data.frame(id=1:length(spol_norm)), FALSE)

# Smaller buffer
spol_small<- gBuffer(sp, width=0.05, byid=TRUE)
spdf_small <- SpatialPolygonsDataFrame(spol_small, data.frame(id=1:length(spol_small)), FALSE)

plot(spdf_norm); par(new=F)
plot(spdf_small)

## Extract values and calculate mean, see results
(ex.mat.norm <- vx$extract(spdf_norm, fun=median))
(ex.mat.small <- vx$extract(spdf_small, fun=median)) # -> 3 NA's

EDIT: Here is the workaround solution by using the raster package for extracting data for the NA polygons and then replacing the NA's with the extracted. When using tons of data, like I do, this solution can speed up the processing quite significantly. Continuing from previous code block:

# Rasterize matrix
r <- raster(mat)
# Use polygon center coordinates (SpatialPoints)
ex.mat.small.temp <- extract(r,sp)
ex.mat.small.temp # check
# Compare
ex.mat.small == ex.mat.small.temp
# Replace NA's with ex.mat.small.temp
require(purrr)
ex.mat.small[,1] <- map2_dbl(ex.mat.small,ex.mat.small.temp, coalesce)

# Fixed outcome
ex.mat.small


texts:
I would like to embed images inside of the toggle controls so that the user can use the toggle menu as a legend, how do I go about this using Mapbox GL JS?

    <!DOCTYPE html>

    <meta charset='utf-8' />
    <title></title>
    <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

    <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.20.0/mapbox-gl.css' rel='stylesheet' />


texts:
When i run the query with ST_Tranform I get this error:

  Input geometry has unknown (0) SRID
  SQL state: XX000

How can I resolve this error? I want to convert from UTM into geographic lat/lon.

This is the query that gives me the error:

select st_asgeojson(st_transform(geom, 4326)) AS geojson FROM ptable where level_1='Built Up'


texts:
OpenlandMap precipitation dataset has a date range between 2007 to 2019. It is an image with the most accurate dataset between precipitation datasets.
Is this possible to filter it's date to a specified date range?
I mean the .filterDate method doesn't work because that's an image, but with a date range.


texts: Is there any gdal tool available for obtaining data from number of files. actually i need  atmospheric parameters which are downloaded from modis website. some datasets are not having data i.e "nodata'. so i have to take previous day's data from database.the files needed are in NETCDf format. i need to check one by one each netcdf file for particular co-ordinate(lat,long) which one is having a value and assign to the variable of the module.can anyone suggest me for obtaining data.


texts: {"Cannot convert value"} error while save geometry data to postgresql table

texts:
I want to add a record to postgresql table with DevArt dotConnect for PostgreSQL component. I use Entity Framework 6 and I add a kind of NuGet packages to my project. SharpMap, NetTopologySuite,GeoAPI...etc common geographic libraries for .net. But i get an error while saving data. Error is {"Cannot convert value"}. I can not solve this problem.

My Code

        using (var db = new MAKSDBDEVEntities())
        {
            var ilce = new MaksIlce()
            {
                Adi = "test",
                Guid = "hede",
                IlkKayitIp = "127.0.0.1",
                IlkKayitTarihi = DateTime.Now,
                IlkKayitKullaniciId = 0,
                UlkeGID = 1,
                Geom = DbGeometry.PointFromText("POINT(10 20)", 4326)
            };
            db.MaksIlces.AddObject(ilce);
            db.SaveChanges();
        }


texts: Why does QGIS Vector OpenStreetMap import lose features?

texts:
I'm attempting to use OSM data in QGIS 2.8.2. I've imported my data into SpatiaLite layers using the following steps:

Vector -> OpenStreetMap -> Download Data...
Vector -> OpenStreetMap -> Import Topology from XML...
Vector -> OpenStreetMap -> Export Topology to SpatiaLite...

For the last step, I imported the polylines. I noticed that some of the primary highways have gaps in the middle. This isn't a case of a segment going outside my region--the gap is right in the middle of the region.

You'll notice that I put the primary highway on its own layer by selecting polylines with attributes I wanted. The light-blue lines beneath the primary highway are every single polyline. Notice that there is no polyline connecting the disjoint highway segments.

I'm assuming that when I chose "Vector -> OpenStreetMap -> Export Topology to SpatiaLite...", all polylines are exported no matter which attributes I check. Is that right? If not, is this a symptom of me not selecting "tunnels", for instance?

This is not a duplicate of this question. The aforementioned question is asking how to import OSM data. I've already demonstrated that I can import the data. The steps I listed are the same as in the accepted answer to 1. My question is about missing polylines.


texts: Clipping raster image in ArcGIS problem column and row

texts:
I am using ArcMap 10.6.

Imagine I have x.tif (100x100 pixels) and Y.tif(10x10 pixels).

If I have to clip x.tif with the extent of Y.tif, the result of Z.tif should be 10x10 pixels. but the result of mine is 10x11 pixels.

The geographic coordinate, datum, and spatial references between x and y is the same. I've also checked the environment on clipping tools and check the "maintain clipping extent".

Is there anything I can do to clip the Z.tif with exact column and row?


texts: Google Earth Engine: How to Refresh a GEE App

texts:
I would like to reset a GEE App to its initial state. That includes removing the added Layers and geometries and clear the user entered values (textboxes,...). The App has a Map and Panel with some widgets. Basically I'm looking for the same behavior as refreshing the web-browser.
I tried to add a ui.Button as follow:

var refreshButton = ui.Button({
  label: 'Reset App!',
  onClick: function() {
  //what to put here??
  }
});

I tried the Map.clear() and Panel.clear()but that is not the behavior I need.
Any ideas ? thanks


texts:
I created a local map service in ArcCatalog.

Then I opened ArcCatalog and in "GIS Servers" clicked in XXX (my computer name). Then listed services down. I right clicked on my service and select services properties. Then I clicked on capabilies. There is a URL TextBox.

Is this my service URL?

If yes: Why This URL not work in browser?

If no: How can I get the service URL?


texts: Applying field mapping during Merge using ArcPy

texts:
I am trying to merge 10 stream feature classes together using arcpy. But I need to limit the number of fields in the output (all feature classes share the same fields).

working_dir = r"C:\Projects\MyProj\stream_lines.gdb"
arcpy.env.workspace = working_dir
all_regions_name = "CNMS_AllRegions_{}.gdb".format(working_dir_name)
all_regions_gdb = os.path.join(working_dir, all_regions_name)

s_studies_for_merge = ['stream1', 'stream2', 'stream3', ... 'stream10']
merged_output = os.path.join(working_dir, "AllStreams")

fields = ['REACH_ID', 'STUDY_ID', 'CID', 'FLD_ZONE', 'VALIDATION_STATUS',
          'STATUS_TYP', 'MILES', 'STATUS_DATE', 'STUDY_TYPE', 'LINE_TYPE',
          'BS_ZONE', 'BS_STDYTYP']
cnms_inventory = os.path.join(all_regions_gdb, "CNMS_Inventory")
s_studs = os.path.join(cnms_inventory, "S_Studies_Ln")

To do this I am trying to apply field mapping during the merge (I specifically need to limit the fields during the merge, not after). I am having difficulty getting only the fields I want mapped and excluding those that do not belong.

fieldMappings = arcpy.FieldMappings()
for item in s_studies_for_merge:
    print(item)
    for field in fields:
        print(field)
        map = arcpy.FieldMap()
        map.addInputField(item, field)

arcpy.Merge_management(s_studies_for_merge, s_studs, mappings)

When I run the merge above, it is still trying to merge a field that I don't want. What am I doing incorrectly here?


texts:
I have added a new SLD style for my nyc_roads layer but no map is displayed with the new style. I get the error "The requested Style can not be used with this layer. The style specifies an attribute of NAME and the layer is: wsNYCRoads:nyc_roads".

The named layer has already been named as suggested but still the same error appears. Why?

 <Name> wsNYCRoads:nyc_roads </Name>


texts: QGIS: Select point(s) in other layer with similar datetime value from selection

texts:
For the interpretation of bird observations, I need to identify the position on a GPS-track (gpx-files), from which the position of a bird observation (other gpx-waypoint layer) was recorded.

The tracks are recorded on an Android device with the Locus app. When a bird is observed, a new point is created on the map in a different point layer.

For the interpretation of the findings later, it would be nice to be able to identify the position from where the observation was taken.

I am looking for ways to develop a tool, where one can select a point in the bird layer, click a button, and the corresponding track points are somehow highlighted.

To start with, I need to find out in general, how I can do a selection query across two layers with the value taken from a selected feature in the first layer.

The next problem would then be, that the timestamps of both layers are never completely the same, as you can see in the example, so in most cases, the exact time of the observation will be between two track points, so a hard comparison of the two layers is not possible.


texts: How to install NetCDF plugin to GeoServer

texts:
From the GeoServer manual:

https://docs.geoserver.org/stable/en/user/extensions/netcdf/netcdf.html

To add a NetCDF data store the user must go to Stores –> Add New Store –> NetCDF.

In my install (Windows, version 2.16, Tomcat), I do not see the NetCDF entry, but only those options:

Raster Data Sources
 ArcGrid - ARC/INFO ASCII GRID Coverage Format
 GeoPackage (mosaic) - GeoPackage mosaic plugin
 GeoTIFF - Tagged Image File Format with Geographic information
 ImageMosaic - Image mosaicking plugin
 WorldImage - A raster file accompanied by a spatial data file

I do not see the way to install the plugin anywhere, and this post Installing netCDF plugin for Geoserver
The link at "Geosolutions IT Training." which may have the plugin installation instruction, is broken.

I also tried installing NetCDF-4 native library https://docs.geoserver.org/stable/en/user/extensions/netcdf-out/nc4.html#nc4 on windows, and it did not help.

There must be some very basics missing here. I do not see the extension/plugin present on my 2.16 GeoServer install, and surprisingly, google around and I do not find anywhere on how to install the plugin.


texts: Moran's I: Selecting the best option for conceptualizing the spatial relationship between # births & county locations

texts:
I'm using Moran's I to see whether there are clusters of births rates present among the 83 Michigan counties, and, if so, whether the clustering is statistically significant.

The inputs are the feature locations (counties) and the feature values (the total fertility rate values we've calculated for each county).

But now I'm stuck because I'm unsure which method I should use to model the spatial relationship. It seems like "polygon continuity", "K nearest neighbors", or the option to "generate spatial weights matrix file" could apply.

Should I select the option to generate a spatial weights matrix, since the total fertility rates for all 83 counties would likely need to be assessed entirely in order to determine the boundaries for clusters?


texts: How do I create a line connecting Features in GeoJSON? I hope to avoid having to re-enter coordinates for the polylines I use to connect the Features separately...


texts: QGIS Batch process of Clipping Raster by Extents using For loop from Python Console

texts:
I am clipping a raster by extents from the python console in QGIS by using the following code :

from processing.tools import *
general.runalg("gdalogr:cliprasterbyextent","E:/Sreeraj/Task/Split/Extract_agri.tif","","-180.0,0.0,0.0,90.0",5,4,75,6,1,False,0,False,"","E:/Sreeraj/Task/Split/output/1.tif")

Here, I am giving extents as minimum_longitude=180.0, minimum_latitude=0.0, maximim_longitude=0.0, and maximum_latitude=90.0.

Now, I want to do a clipping using processing.runalg in QGIS python console, where I will use a for loop for clipping 4 output .tif files as a batch process.

For doing this batch process of clipping, I tried the given below code:

output_location = "E:/Sreeraj/Task/Split/output/"
output file = 1
for long in range(-180,0,5):
  longmax = long + 5
  for lat in range(0,90,5):
    latmax = lat + 5
    extents = "long,lat,longmax,latmax"
    output = str(output_location) + str(output file) + ".tif"
    general.runalg("gdalogr:cliprasterbyextent","E:/Sreeraj/Task Global Data CONVERSION/Split/Extract_agri.TIF","",extents,5,4,75,6,1,False,0,False,"",output)
    output file = output file + 1

Unfortunately, I am getting error message as :

"Error: Wrong parameter value: extents"

Instead of directly giving the extent values like "-180.0,0.0,0.0,90.0" inside general.runalg(), how can I give  "extents" (if extents = "-180,0,0,90").

I think that the issue is with :

extents = "long,lat,longmax,latmax"

When I do like this :

extents = "-180.0,0.0,0.0,90.0"
general.runalg("gdalogr:cliprasterbyextent","E:/Sreeraj/Task Global Data CONVERSION/Split/Extract_agri.TIF","",extents,5,4,75,6,1,False,0,False,"",output)

then the code is working.

But, if I do like this :

long = -180.0
lat = 0.0
longmax = 0.0
latmax = 90.0
extents = "long,lat,longmax,latmax"
general.runalg("gdalogr:cliprasterbyextent","E:/Sreeraj/Task Global Data CONVERSION/Split/Extract_agri.TIF","",extents,5,4,75,6,1,False,0,False,"",output)

then I am getting the error.

So, clearly, the issue is with :

extents = "longmin,longmax,latmin,latmax"

Please help me to solve this error.


texts: QGIS: static Atlas, it do not follow the point-layer

texts:
I do not use Atlas for years.
Now I set Atlas to follow a point layer but when I export in separate files it gives me this error:
impossible write on C:...path\name_item.pdf because the file could be used by another application or path not valid (I try with different file names and different folders)
When I try to export it as a single file it works
I start with a new project but te problem was not solved


texts: Snapping point to point nearest to it in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.3 for Desktop.

I have 2 point features centroid and station. I want to snap all the points of centroid to the nearest point in station such that centroid is only snapped with points of station.


texts: How to use r.mapcalculator to change values within a raster cost map to negative

texts: I have created a cost map in r.cost module in GRASS (through QGIS). This is actually a distance map from a protected area and distance from protected area has been calculated in both directions. How I can change the values within protected area to negative? I think I should use r.mapcalculator module in GRASS.


texts:
Has anyone seen the following problem?

I've created a sample ArcGlobe application with Visual Studio - this obviously gives me the template where I should be able to launch the app.  Said application launches, I then go to load a .3DD File from the preexisting samples that come with ArcGIS.  I try loading the .3DD file, the APP hangs for a bit and eventually the Globe Control comes up.

When this happens, I see the World Physical Map (layer ) and the US Elevation (layer).  But the weird thing is is that the screen is black.  If I uncheck the layers, then I get a wireframe globe with no data displayed.

Any thoughts?


texts: Could not figure out OD cost matrix message

texts:
I have street nodes as the origin points (I have 19,133 for the whole area) and grocery stores(I have 9 of them in all) as the destination points.  After clicking on the solve button for the OD cost matrix analysis, I get a message display as shown in the sample that I attached.  (Not all origins get connected by straight lines to the destination, I get only a few of these origins that connect)

Anyone could help me out with this issue?

@fatih_dur I tried using the topology toolset to correct errors in the road network.  After running the error inspector, I found out that there are 116 errors.  I was trying to correct using the "Fix Error Tool" but I could not proceed.  Could you help me out?


texts:
how to convert a map (.sqlite) to RMaps by MOBAC or other app?

THE MAP TYPE : OSMDroid
MAP FILE NAME : NEW.sqlite
MAP ATTACHMENT FILES : NEW.sqlite-journal

I don't have more experience about this
so please i need easiest way to do this with simple explanation


texts: Adding feature to PostGIS layer that has joins is really slow

texts:
I have a large project that I've been setting up to work with QGIS 3.12 and I have now hit a problem with using "joins" to join several attribute tables with a geometric one.

I have reduced the problem to a minimal setup that has only two layers with a one to one (or zero) relation:

point_layer - a layer with a geometry (points) that has 64k entries
attribute_layer - a layer with one text field and a foreign key to the point_layer

On the point_layer I added a join to the attribute_layer so that we can edit the attributes directly in the feature form.

After doing this whenever I try to create a point on the point_layer it takes more than 1 min to create the feature.

By examining the PostgreSQL log I found out that QGIS is querying each feature on the point_layer for it's relation on the attribute_layer. These queries are made one by one, so even though I'm running the database on my own computer, the sheer number of queries leads to very poor performance.

I've already tried tweaking the cache config on the joins, hiding the form on add feature and enabling the "automatically create transaction groups where possible" to no avail.

I've also tested older QGIS versions 3.4 and 3.8 without success.

Am I doing anything wrong? Is there any way to work around this?

UPDATE:
As requested by vince here is the table schema:

CREATE TABLE point_layer (
    ogc_fid integer PRIMARY KEY,
    geom public.geometry(GeometryZ,3763),
    height double precision
);

CREATE INDEX point_layer_geom_idx ON point_layer USING gist (geom);

CREATE TABLE attribute_layer (
    id SERIAL PRIMARY KEY,
    fid integer REFERENCES point_layer(ogc_fid),
    some_value text
);

The query that is made repeatedly is (the fid changes):

SELECT "id","some_value"::text,"fid"::text
FROM "public"."attribute_layer"
WHERE ("fid" = ?) LIMIT 1

The explain analyze is (the table is empty):

                                                  QUERY PLAN
---------------------------------------------------------------------------------------------------------------
 Limit  (cost=0.00..4.17 rows=1 width=68) (actual time=0.014..0.014 rows=0 loops=1)
   ->  Seq Scan on test_relation  (cost=0.00..25.03 rows=6 width=68) (actual time=0.011..0.011 rows=0 loops=1)
         Filter: (fid = 29853)
 Planning Time: 0.506 ms
 Execution Time: 0.086 ms


texts: Where can I programmatically download SRTM_V4.1 data?

texts:
SRTM v4 is an effort by the folks at CIAT-CSI to scrub, polish, and remove gaps in the raw SRTM data released by NASA (SRTM 3).

One userfriendly interface : http://dwtkns.com/srtm/ provides srtm_xx_yy.zip tiles of ~20Mo each.

The CIAT-CSI project provides 3 zipped tiles (W, NE, SE) of ~1Go each. These files are behind a public password (see page), on box.com, which block programmatic download.

Is there any mirror with these 3 big SRTM 4 files using a stable link ? so I may use a makefile to let coworkers to programmatically download them via curl or alikes.


texts: Copying fields between two layers in ArcGIS Desktop?

texts:
I have this attribute table (part of the whole attribute table) from one layer:

and this attribute table (part of the whole attribute table) from the second layer:

I try copy the cells value in fields: name2, source into the second layer:

Is there a way to do it in ArcGIS Desktop 10.3?

i can't do join because those two layers have no field with the same values.
i also red

Copying features from one layer into different layer using ArcGIS for Desktop?
Copy table to attribute table


texts:
First, I created a regular grid of the size 1km^2 (epsg: 2056). Then I calculated the area of of this grid using different GIS:

Python: arcpy

  area_arc

# Python Setup
import arcpy
arcpy.env.overwriteOutput = True

# Data_in
poly = "C:/Hexbins.shp"

# Area Caculation
arcpy.AddField_management(poly,'area_arc','DOUBLE')
arcpy.CalculateField_management(poly,'area_arc','!shape.area!','PYTHON')
#arcpy.CalculateGeometryAttributes_management(poly,["area_arcpy", "AREA"],'Meters')

R sf

  area_sf

Hexbins<- read_sf("LC:/Hexbins.shp")

Hexbins<-Hexbins %>%
mutate(
area_sf= as.vector(st_area(.))
)

QGIS

  area_q

Field calculator:
Input: poly
- output field name: area_q
Output field type: Decimal number (real)
Outpu field length: 10
Precision: 3
Expression: $area

  area_q2

Field calculator:
Input: poly
- output field name: area_q
Output field type: Whole number
Outpu field length: 10
Precision: 3
Expression: $area

The libraries Python arcpy and R sf lead to the expected result of 1km^2. The result I obtained using qgis is somewhat astonishing:

The libraries Python arcpy and R sf lead to the expected result of 1km^2. The result I obtained using QGIS $area is somewhat astonishing:

Is there an explanation for this randomness?


texts: Syncing north arrow with polygon centerpoint

texts:
I am using the QGIS 3.2.1 atlas functionality to create a series of maps. The maps are of polygons from 60 degrees north to 78 degrees north. The north arrow seem to be aligned based on its position within the map.

I am looking for a way to connect the arrov direction with the respective polygon presented on the actual atlas plate. I would probably need a fuction which creates a value for the offset based on a calculated longitudal senterpoint.

Is this possible? If so - is there code available?


texts: Geotools & Natural Earth finding lat/lon not always accurate

texts:
I am trying to take a lat/lon and find which country the lat/lon is located in.

Below is a boilerplate of code I found to accomplish this, which was found in this post (Reverse Geocode latitude longitude with Geotools using Natural Earth Shapefiles in Java
)

SpatialIndexFeatureCollection countries;
final static FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
final static public String countryFile ="/data/natural_earth/10m_cultural/10m_admin_0_countries.shp";

public SimpleGeoCoder() throws IOException {
  // load the country shapefile
  URL countryURL = DataUtilities.fileToURL(new File(countryFile));
  HashMap<String, Object> params = new HashMap<>();
  params.put("url", countryURL);
  DataStore ds = DataStoreFinder.getDataStore(params);
  if (ds == null) {
    throw new IOException("couldn't open " + params.get("url"));
  }
  Name name = ds.getNames().get(0);
  countries = new SpatialIndexFeatureCollection(ds.getFeatureSource(name).getFeatures());
}

GeometryFactory gf = new GeometryFactory();
Point london = gf.createPoint(new Coordinate(0.0, 51.0));

public SimpleFeatureCollection lookup(Point p) {
  Filter f = ff.contains(ff.property("the_geom"), ff.literal(p));
  return countries.subCollection(f);
 }

SimpleFeatureIterator itr = features.features();
try {
  while (itr.hasNext()) {
    SimpleFeature f = itr.next();
    System.out.println(f.getAttribute("NAME"));
  }
} finally {
  itr.close();
}

Most of the time this solution works, however there seems to be cases where a lat/lon will not be found within the Natural Earth shapefiles.  A perfect example would be for the city of Miami Florida (25.761681, -80.191788) when using the shapefile "10m_cultural/10m_admin_0_countries.shp".

My original solution was to use the shapefile "110m_cultural/110m_admin_0_countries.shp" the Miami example will be fine.  However I have found a few other examples which occur within the 110m shapefile and I don't know what the solution is.  Do I use both the 10m and the 110m shapefiles?  Or is there some setting or something else I need to do to accurately solve this issue?


texts: Creating an Action in QGIS with multiple Expression?

texts:
I want to create an action in QGIS with multiple Expression. But unable to do so. I want to open images from two different directories with one action. Is it possible? How can I do so?

e.g.

Type : open
name : open Photo
D:\image\Path_1\[%"Photo"%].jpg
D:\image\Path_2\[%"Photo"%].jpg

But it's not working. How can I fix it?


texts: Finding city polygon for Brisbane, Australia using OpenStreetMap?

texts:
I can't seem to find a polygon for the city of Brisbane, Australia. I can only find a Point. Is there another name i should be searching for? ie. city, municipality.

I've searched for "Brisbane, Australia", "Brisbane City, Australia", "Brisbane Municipality, Australia", "Brisbane Region, Australia"... I'm looking for a polygon that circles the city of Brisbane. Similar to what comes up when you search for "Greater Melbourne, Australia" (place_id=156135873) or "Sydney, Australia" (place_id=156208530). It needs to come from OpenStreetMap.


texts: How do I figure out the length of polyline features in a feature class with a WGS 1984 coordinate system?

texts:
I'm trying to figure out lengths of polyline features in a feature class with a WGS 1984 coordinate system.

Here's code to list lengths:

>>> import arcpy
>>> inFc = r"C:\Users\e1b8\Desktop\E1B8\GIS_Stackexchange\data.gdb\test"
>>> sr = arcpy.Describe (inFc).spatialReference
>>> sr.name
u'GCS_WGS_1984'
>>> with arcpy.da.SearchCursor (inFc, "SHAPE@") as cursor:
    for geom, in cursor:
        print geom.length

0.000785568606752
0.000784558202405
0.00117166138222
0.000816779544786

I'm not sure what units these lengths are in. I know it's generally in the coordinate system's linear units, but WGS 1984 has no linear units.

>>> sr.linearUnitName
u''

Ultimately I'd like to be able to write code that works with distances and is compatible with both projected and geographic coordinate systems. Usually I use the spatial reference's metersPerUnit property to get to a common unit regardless of the coordinate system's unit, but this throws an error for WGS 1984:

>>> sr.metersPerUnit

Traceback (most recent call last):
  File "<pyshell#13>", line 1, in <module>
    sr.metersPerUnit
AttributeError: SpatialReference: Get attribute metersPerUnit does not exist

So again, how can I get the length of features? And preferably, how do I figure out their lengths in meters?


texts:
How might I go about exporting a polygon shapefile (e.g. CSV) so that attributes (not geometry) are editable by others, and importing it with the changed attributes (of which there are several columns worth)?

I tried the mmqgis plugin (Transfer > Attributes Export/Join to/from CSV file), but
1) it seems to only be able to import one attribute at a time and
2) creates a separate temp.shp without geometry and
3) does not attach the new attributes to the existing polygon shapefile.

Is there something I'm doing wrong here, or do you have a better solution for my task at hand?

Thanks!