texts: Get the attribute value of another layer with geometry condition

texts:
With QGIS 3.14 and the field calculator :
I have two layers, one with stream network and one with watersheds.
I would like to attribute the name of the watershed to the streams contained inside it with an expression in the field calculator. How to include the "within" condition ?
I guess it should be something like that
CASE
WHEN  within(  $geometry , $geometry) THEN
get_feature('watersheds',"name", $currentfeature )
END

But replacing the second $geometry with a reference to the polygons of the watershed layer
With PyQGIS 3
I also tried the following PyQGIS code , that requires some more informations.
streamInput
The layer containing the stream network has one unique ID per each stream ('IDjoin2')
The layer of watersheds contains several "levels" of watersheds. The main streams belongs to the bigger watersheds, and the smallest have their own little watershed, but also belongs to bigger ones. There are 7 levels of watersheds : the 1st one is the bigger one, the last one is the smallest.
streamBvInput
A spatial join (1-to-many) with the watersheds layer allowed me to have, for each stream, all the watersheds in which it is contained.
csv_path
I need to get, for each stream, the name of the smallest watershed that contains it, that's why I calculated statistics on the above layer, to obtain the maximum level of watershed for each IDjoin2/stream.
project = QgsProject.instance()

#input layers
streamInput = QgsProject.instance().mapLayersByName('Reseau_Hydro_..._Nom')[0]
streamBvInput = QgsProject.instance().mapLayersByName('Reseau_Hydro_..._Nom_intersectBVmulti')[0]
csv_path = 'C:/Users/.../IDjoin2_intesect_maxBV.csv'

streamLayer = project.mapLayer(streamInput.id())
streamBvLayer = project.mapLayer(streamBvInput.id())

if not (streamLayer.isEditable()):
    streamLayer.startEditing()

layerStreamFields = streamLayer.fields()
with open(csv_path, newline='') as csvfile:
    reader = csv.DictReader(csvfile)
    for row in reader: #broswe into statistics csv file
    # get the highest level of watershed for each stream
        IDcsv = int(row['IDjoin2'])
        searchBVmax = row['max']

        iterStreamBV = streamBvLayer.getFeatures()
        for featureBV in iterStreamBV: #broswe into the join file
            searchIDbv = featureBV['IDjoin2']
            if (IDcsv == searchIDbv): # if ID in csv and bvInput are identical
                if (featureBV['NIVEAU_BAS'] == searchBVmax) :
                #if watershed level in csv and bvInput are identical
                    BVname = featureBV['NOM_COURS_1']
                    # get the name of the watershed
                    iterStreamFeatures = streamLayer.getFeatures()
                    for featureS in iterStreamFeatures:
                        searchIDj = featureS['IDjoin2']
                        if (searchIDj == searchIDbv):
                            featureS.setAttribute('NomBV2',BVname)
                            attribute the name to the stream in streamInput

                        streamLayer.updateFeature(featureParcelle)

streamLayer.commitChanges()

But as soon as I run this script QGIS is freezing, and I received no error messages. I began testing the code almost line per line, and it looks like a part of the problem comes from the size of the layer "streamBvInput", that contains more than 114 000 features.


texts: Applying a definition query to date field for feature in file geodatabase

texts:
Is it possible to apply a definition query to show data for today, yesterday and two days ago for a field where the data type is 'Date'?

I have a script that downloads MODIS fire data and wish to show data on a map for only the past three days, but am totally stumped on how to write the definition query. It should be something like SELECT * FROM MODIS WHERE: DATE <= CURRENT_DATE, but how to get the query down to only three days has got me stumped.

It would be nice to not have to manually type the values of the date in every time someone asks where a fire is burning near town. I also realize the MODIS data is available which shows only the past seven days, but showing fires from only three days is the goal.

The date data is formatted as date '2015-01-12 00:00:00' in the query builder.


texts: ArcObjects ITable.Search hits default table version but cursor iterates over edited current version of table?

texts:
I'm very new to ArcObjects and working on an ArcMap addin. I'm having some trouble with understanding the results I'm receiving from my queryFilter. It seems doing a ITable.Search(queryFilter, true) [recycling set to false nets the same results] hits the default table version but returns a cursor that iterates over the current version (active in ArcMap) of the table.

In ArcMap when I create a new version of my geodatabase and run my addin tool I get correct results. However, if I make an edit in the centerline_sequence table in ArcMap and then run the tool again. It appears to run the query against the default or parent version as it returns rows that should have been excluded by my edit. When the cursor cycles through the results I get the edited rows which should have been excluded from the results due to the query select statement. The real oddity to me is the cursor rows include the record edits which would have excluded that record from the results.

I have referenced the table originally like this:

IWorkspace workspace = ((IDataset)map.Layer[0]).Workspace;
featureWorkspace = (IFeatureWorkspace)workspace;
centerlineSequenceTable = featureWorkspace.OpenTable(tableName);

Then thinking I wasn't getting the current version I got the table this way:

IFeatureWorkspace cws = (IFeatureWorkspace)this.childVersion;
centerlineSequenceTable = cws.OpenTable(tableName);

However both showed they were the same active version in ArcMap.

In my queryFilter.WhereClause I pass a statement that looks for a records which include a specific field value. I have tested this clause in SQL Developer and believe it is returning correct results.

I search the centerlineSequenceTable like:

IQueryFilter2 queryFilter = new QueryFilterClass();
queryFilter.WhereClause = whereString;
queryFilter.SubFields = objectIdField + "," + routeIdField + "," + roadwayIdGuidField + "," + dominanceValueField;
IQueryFilterDefinition queryFilterDefinition = (IQueryFilterDefinition)queryFilter;
queryFilterDefinition.PostfixClause = "ORDER BY " + roadwayIdGuidField + ", " + routeIdField;
ICursor cursor = centerlineSequenceTable.Search(queryFilter, true);

Gather my results like so:

IRow row = cursor.NextRow();
while (row != null)
{
…
row = cursor.NextRow();
}

I'm not much of an ArcGIS user so quite possibly I'm just expecting something that isn't going to happen in my workflow without reconciling the current version or I'm missing something in how I reference and/or search the table.


texts: Project raster when some of extent out of range

texts:
I have a raster that is in a North America Albers equal area projection. I would like to project this raster to unprojected long/lat, but I get an error. I understand that this is likely because some of the northern-most cells would transform to > 90 degrees, but I haven't been able to think of a way to work around this problem.

I've put together a reproducible example:

library(raster)
library(maptools)
data(wrld_simpl)

# generate a raster of north america
bb <- rbind(c(-166, 78), c(-175, 48), c(-108, 0), c(-55, 12), c(5, 78), c(-166, 78))
bb <- SpatialPoints(bb, proj4string = CRS('+proj=longlat +datum=WGS84'))
bb <- rgeos::gConvexHull(bb)

rasterTemplate <- raster(ext = extent(wrld_simpl), res = c(1,1))
rasterTemplate <- rasterize(wrld_simpl, rasterTemplate)
rasterTemplate[!is.na(rasterTemplate)] <- 1

rasterTemplate <- crop(rasterTemplate, bb)
rasterTemplate <- mask(rasterTemplate, bb)
rasterTemplate <- trim(rasterTemplate)

# here is our equal area projection
EAproj <- '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

# project to the EA projection
e <- projectExtent(rasterTemplate, crs = EAproj)
rasterTemplate <- projectRaster(rasterTemplate, to = e)
rasterTemplate <- trim(rasterTemplate)

plot(rasterTemplate, legend = FALSE)

# project raster to long/lat
e <- projectExtent(rasterTemplate, crs = '+proj=longlat +datum=WGS84')

Error in if (maxy == miny) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
  36 projected point(s) not finite

rasterTemplateLongLat <- projectRaster(rasterTemplate, to = e)

Any suggestions?


texts: Prevent the Feature Attributes dialog box in QGIS

texts: Any time I am inserting new point feature in QGIS, the Feature Attribute dialog box is popping up. Is there a way to prevent it? (I am willing to insert many points at once and this pop up box slows down the process).


texts:
I am working on NDVI in ENVI... the output file is in Envi Standard format... I required it in Envi Classification format, which can be done through Envi Edit Header tool... This requirement is for its use in Confusion Matrix using Ground Truth Image... but when i performed this, the resulted NDVI image transformed into BLACK color showing as 'unclassified OR class1' portion (what ever name was suggested at the time of conversion into classification format). (As by the editing header, from standard format to classification, it require at-least 2 classes). I need this to perform Tasseled Cap comparison with NDVI.

I was not given the option to change the format of output... either integer or other way... It only give this error when i tried to change the output format from Envi Standard to Envi Classification by Edit Header tool.

This error is due to heavy file size... most of the time as a glitch in ENVI outputs.
Although the process runs smoothly with "ENVI Edit Header" tool... it automatically changes the output into Classification format with as many classess as defined by the user.


texts: Is there a way to install leaflet search plugin without using npm?

texts:
I am not familiar with Node.js or npm at all, but my IT has told me node.js/npm is not compatible with our IIS, and the only search plugins to add on to leaflet are using npm to install them.

I did a search here and found this solution without npm - leaflet-control-geocoder, but it searches address databases such as google or bing, not the geojson data I have plugged into my map that I need to be able to search.

So my question is, is there a way to use the Leaflet Search plugin or leaflet-fusesearch plugin without using npm to install it?


texts: How to set srid of a single point coordinates pair?

texts:
I think is better explain with a small example:
selected a geometry point i get it back in a coordinate pair

SELECT ST_AsGeoJSON('010100002031BF0D009A999919F9BD2D41CDCCCC4C758E5241')
--result:
"{"type":"Point","coordinates":[974588.55,4864469.2]}"

to return back the geometry simply i make

SELECT ST_GeometryFromText('POINT(974588.55 4864469.2)')

having for example the same pair of coordinates but in a different "form" (srid?)

coordinates:[9.7458855,48.644692]}

How i can transform them into the previous pair?
I know,i'm a bit confused about all this stuff, where i can find a good reading source?


texts: shp2pgsql uses deprecated addGeometryColumn function

texts:
It is very likely that I am misunderstanding the documentation as I'm new to PostGIS, but it says here regarding the "addGeometryColumn" function that since version 2.0.0 of PostGIS:

  This function no longer updates geometry_columns [...] So for example building a wgs84 POINT column with this function is now equivalent to: ALTER TABLE some_table ADD COLUMN geom geometry(Point,4326);

However, the command line tool "shp2pgsql" which was shipped with the "postgis" package (I am on Ubuntu 16.04 LTS xenial) generates SQL containing the addGeometryColumn() function, for example:

SELECT AddGeometryColumn('','mytestgeo','geom','26918','POINT',2);

I am definitely using PostGIS 2.5:

SELECT postgis_version();
2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

Since the shp2pgsql tool was shipped with the postgis package, I would have expected it to conform to PostGIS 2.5 specification. Why then does it still use this method? Is it for backwards compatibility? Have I misunderstood the docs?


texts:
I have a small table of 10 cities (in one table) that I am snapping to the nearest road (in another table) with the following select and storing the table snapped_cities:

Select b.id, ST_Closestpoint(ST_Collect(a.geom), b.geom) as geom
into snapped_cities
from roads a inner join cities b
on ST_Dwithin(a.geom, b.geom, 1000)
group by b.id;

I would however like to have road ids and not just the geom (a.id above). How do I modify the select above to get these?

I use postgresql 9.2 and postgis 2.0

thanks,


texts: Add scaleline control to OpenLayers 3 map with unknown projection

texts:
I've a MapGuide map display with OpenLayers 3 like this:

var options = {
"minZoom": 0,
"maxZoom": 100,
"center": [mapInfos.center.x, mapInfos.center.y],
"extent": [mapInfos.extent.minX, mapInfos.extent.minY,   mapInfos.extent.maxX, mapInfos.extent.maxY]
};
var projectionUnit = mapInfos.projection.unit;
var projectionCode = mapInfos.projection.code;
var projection = new ol.proj.Projection({
  code: projectionCode,
  extent: [mapInfos.projection.extent.minX, mapInfos.projection.extent.minY, mapInfos.projection.extent.maxX, mapInfos.projection.extent.maxY],
  units: projectionUnit,
});
ol.proj.addProjection(projection);
options.projection = projectionCode
var imgMG = new ol.source.ImageMapGuide({
  url: me._agentUrl,
  useOverlay: true,
  useAsyncOverlay: true,
  metersPerUnit: me._metersPerUnit,
  params: {
    FORMAT: 'PNG',
    MAPNAME: me._mapName,
    behavior: 2,
    SESSION: me._session
  },
  ratio: 2
});
var layer = new ol.layer.Image({
  source: imgMG,
});
me._map = new ol.Map({
  interactions: ol.interaction.defaults({
    doubleClickZoom: false,
    shiftDragZoom: false,
    altShiftDragRotate: false,
    keyboard: false,
    pinchRotate: false,
    pinchZoom: false,
  }),
  controls: [],
  layers: [layer],
  target: me._options.target,
  view: new ol.View(options)
});

And when I want to add a scaleline control to my map like this:

var scaleLineControl = new ol.control.ScaleLine({
  "className": this._control.className,
  "units": "metric" //To do dynamically with the map projection units /!\
});

I've an error: "transform should be defined"

Have you an idea to resolve my problem?


texts: GDAL (Python) Aggregate Raster Into Lower Resolution

texts:
I have a (continent-scale) raster with 30m resolution in an Albers equal area projection. Each 30m cell has a value 1-34.

I would like to aggregate information from this raster into a 1km grid such that each 1km grid cell contains a floating-point value which is the number of 30m cells with a value of 2 divided by the total number of 30m cells which fit into the 1km cell.

Is there a way of doing this with GDAL? Or an efficient way to hook GDAL and Python together to perform the calculation (given the large size of the raster, not all of it fits into RAM at once)?


texts: Standalone PyQGIS script for QGIS 3.8 processing

texts:
Using the code examples for the processing module found in the cheatsheet, in the QGIS Python console I can print a list of the algorithms available. I want to do the same from a "standalone" python script for my Windows box. Is it possible?

I made a Windows batch script based on "C:\OSGeo4W64\bin\python-qgis.bat" (from this answer) to setup the environment:

@ECHO OFF

REM this batch file sets required environment variables to allow pyqgis scripts to work
REM and then runs VS Code in that environment.

set OSGEO4W_ROOT=C:\OSGeo4W64

@echo off
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"
call "%OSGEO4W_ROOT%\bin\qt5_env.bat"
call "%OSGEO4W_ROOT%\bin\py3_env.bat"

@echo off
path %OSGEO4W_ROOT%\bin;%OSGEO4W_ROOT%\apps\qgis\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis
set GDAL_FILENAME_IS_UTF8=YES

set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis\python;%PYTHONPATH%

cd /d %~dp0

SET VSCBIN=C:\Program Files\Microsoft VS Code\

call "%VSCBIN%\Code.exe" %*

Then after I run that batch script I create a python script in VS Code with the same content as in the cheatsheet:

from qgis.core import QgsApplication

for alg in QgsApplication.processingRegistry().algorithms():
    print("{}:{} --> {}".format(alg.provider().name(), alg.name(), alg.displayName()))

But it does not produce any output. It doesn't produce any errors so that's one small bit of progress. This code also doesn't work when run from the python console opened from running the batch file "C:\OSGeo4W64\bin\python-qgis.bat".

I can see from other examples here and there that I might also need to init QgsApplication like so:

from qgis.core import QgsApplication

# Supply path to qgis install location
QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis", True)
# (not sure why the above line is needed as the batch file already sets the prefix path)

# init Qgis, but don't create the GUI
qgs = QgsApplication([], False)
qgs.initQgis()

for alg in QgsApplication.processingRegistry().algorithms():
    print("{}:{} --> {}".format(alg.provider().name(), alg.name(), alg.displayName()))

But still no output. I get no errors, so I know that the imports must be fine. There just doesn't seem to be any "link" to QGIS itself. What am I doing wrong?


texts: What is MHHW Mean Highest High Water, in terms of datums/coordinate systems?

texts:
I'm modeling sea level rise, which entails adding the projected sea level rise for a specific decade to the current MHHW height (I think. I am new to this.)

My question is very general--what is MHHW in terms of datums? for example, WGS84?

I have projected Sea level rise in centimeters, and have created a raster from it. But I need to add it to the MHHW for the region.

Where do I get MHHW? Do I look at tide charts and average them? The procedure I'm following seems to indicate that MHHW should be considered 0, or just a baseline, ground level.

How can this be right? isn't MHHW above 0? Also I know in some systems sea level is not defined as 0.


texts: How can I work out a percentage cover of a polygon layer by another?

texts:
I am working with ArcGIS pro 2.2.0.

I'm trying to work out how much area is affected by some ocean exclusion areas, I have a grid of 10x10km squares and wanted to work out how much each square is affected by these exclusions.

This is my overall grid area: (red square is just showing location of the second image)

This is an example of the exclusions I want to calculate how much of each grid cell is affected by the exclusions.


texts: Images in my Click Info box are not appearing in both published CartoDB layers

texts:
I'm using the simple web version of CartoDB. I uploaded two datasets, each is a layer and both have the same image URL type (URLS for .png Twitter account images). When I visualise the maps individually, the images appear fine so I know the problem is not with the URLs themselves, but for reasons unknown only one set of images (blue points) displays in the published map - see here: https://drarnakarick.cartodb.com/viz/de49f3ca-cbac-11e5-a9f3-0ecd1babdde5/embed_map

How do I make it so both layers display properly. I've chosen layers because the content/columns differ between datasets.

The custom HTML for the info window (format is the same for both):

x

Company Name

{{compname}}

Headquarters

{{headquarters}}

Description

{{description}}

Official Website

{{website}}

Twitter Profile

{{twitterhandle}}

Other Locations

{{otherlocations}}


texts: How do I access data for each pixel within an image and store it in an array in Google Earth Engine?

texts:
I want to create an array in Earth Engine that is composed of the NDVI data of each pixel in the image. Below is the code that calculates the NDVI of each pixel in the ROI (a region). I can store the data for a specified pixel (check the code below the "Extract the data" comment) in a variable(dataN), but I am unable to store the data for each pixel in the region in an array.

var L8 = ee.ImageCollection("LANDSAT/LC8_L1T_TOA"),
roi = /* color: #d63000 */ee.Geometry.Point([74.47821320074036, 31.556778358715253]);

var rgb_vis = {min: 0, max: 0.3, bands: ['B4', 'B3', 'B2']};
var filtered = L8.filterDate('2015-09-01', '2015-10-30').filterBounds(roi);
var image = ee.Image(filtered.first());
var red = image.select('B4');var nir = image.select('B5');

var ndvi = nir.subtract(red).divide(nir.add(red)); //layer from which I want to access the individual pixel data from

Map.addLayer(image, rgb_vis, 'RGB');
Map.addLayer(ndvi, {min: 0, max: 1}, 'NDVI');

// Extract the data
var data = ndvi
.select("B5")
.reduceRegion(ee.Reducer.first(), roi, 10)
.get("B5")

// Convert to Number for further use
var dataN = ee.Number(data)

// Show data
print(dataN);


texts: Knowing if there is correlation between roads and observations in ArcMap?

texts:
I am using ArcMap 10.3.1.

I am trying to find out whether some observations have been made randomly throughout a landscape, or they are biased towards the location of roads and paths. I have a layer with the observations and a layer with the roads.

Any tip on how I may do this, other than through visual inspection?


texts: Computing Shannon-Weiner biodiversity index using QGIS?

texts:
I wonder if anybody has some experience in calculating biodiversity index (such as Shannon-Weiner) using QGIS?

I am also interested for similar experiences using R scripts within QGIS.

The analysis will use species count from individual stations of marine scientific surveys.

I have in mind to calculate the Shannon index of species diversity for each station of the sampling survey and then to "allocate" the average index within each cell of a vector grid or by interpolation in the Area of interest.

I wonder if anyone has experience in doing this using a spatial query (providing  stations and vector grid are stored in a sqlite database)?


texts: Inserting geometry into spatialite - no geometry

texts:
I have created a test spatialite database and I am having trouble inserting geometry into the table. I have successfully created the geometry column in the table and then I run the insert command. The weird thing is that it runs without error, and the attributes get passed in, however the geometry column is empty when I view it through a ODBC connection, as well as no geometry coming through in QGIS.

My code to create the table:

with sqlite3.connect(r'somedatabase.sqlite') as conn:

    conn.enable_load_extension(True)

    conn.execute("SELECT load_extension('mod_spatialite')")
    conn.execute("SELECT InitSpatialMetaData()")

    c = conn.cursor()

    c.execute("""CREATE TABLE IF NOT EXISTS gid(gid_id INTEGER NOT NULL PRIMARY KEY AUTOINCREMENT,
                latitude REAL, longitude REAL)""")
    c.execute("SELECT AddGeometryColumn('gid','the_geom', 4283,'POINT','XY')")
    conn.commit()

And this is the code to create the geometry

someGeometry = [(150.895,-33.961),(150.896,-33.959),(150.893,-33.960)]
for longitude,latitude in someGeometry:
    c.execute("""INSERT INTO gid(gid_id,latitude,longitude,the_geom)
                 VALUES(NULL,?,?,GeomFromText('POINT((%f %f)',4283))""" %(longitude,latitude),
                    (latitude,longitude,))

conn.commit()

I have also tried to remove the python string formatting by inserting only 1 point with specific coordinates like this

c.execute("""INSERT INTO gid(gid_id,latitude,longitude,the_geom)
                 VALUES(NULL,?,?,GeomFromText('POINT((150.895,-33.961)',4283))""",(latitude,longitude,))


texts: SIGABRT when calling QgsMeshLayer class (QGIS)

texts:
I am working on a C++ application that uses the QGIS API to overlay weather data from NetCDF files onto a KDE Marble map.  I am using version 5.11.3 of Qt and 3.8 of QGIS on Fedora 28.  When trying to load the data using the application, I receive an abort signal (core dumped) and the following error:
pure virtual method called
terminate called without an active exception
__GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:50
50        return ret;

The function call that triggers this in my program is the following:
QgsMeshLayer *vl_main = new QgsMeshLayer("path/to/file.nc","met_data","mdal");

yet I call it just as the qgsmeshlayer.h file states in comments:
 * \subsection mdal MDAL data provider (mdal)
 *
 * Accesses data using the MDAL drivers (https://github.com/lutraconsulting/MDAL). The url
 * is the MDAL connection string. QGIS must be built with MDAL support to allow this provider.

 * \code
 *     QString uri = "test/land.2dm";
 *     QgsMeshLayer *scratchLayer = new QgsMeshLayer(uri, "My Scratch Layer",  "mdal");
 * \endcode

I have ensured that QGIS was built with MDAL support and calling this function is the same method that QgsVector and QgsRaster layers are added, which do not return the pure virtual message.
$ grep "MDAL" CMakeCache.txt
//Determines whether MDAL support should be built
WITH_INTERNAL_MDAL:BOOL=TRUE

To add to the confusion, I haved initialized QGIS in the main application program main.cpp with a QgsApplication object and have all the QGIS plugins and dependencies added to path.  Further, I can add the mesh data when using QGIS instead of QT Creator; however the project needs to work in QT Creator with the other applications.

If any other information to help solve the root of the SIGABRT message is needed from my end, I can append this entry further.  I should note that this version of QGIS was built from source (github) and as it is newer than what is supported through the EOL Fedora 28 RPMs, required newer dependencies to be built, so it could be a version compatibility issue as library conflicts have risen in the past.


texts:
I’m looking for a way to easily trace a transmission line for a certain distance.
 Here is the scenario that I am working with. We have an issue with our transmission line; our OMS tells us that there is a problem “X” number of miles from substation “Y”.
I want to be able to start from that substation and trace that transmission line upstream for the number of miles indicated to the location of the issue so we can dispatch a crew to that location.

I was hoping that someone had something that worked similar to the “Utility Network Analyst” tool.

Any ideas?


texts: How to plot a common colorbar for mutiple raster plots in a single image in R?

texts:
I am plotting a spatial map of monthly temperature. I need to keep a common scale for all the months. How to do this in R?

I have included my code below with the generated map. How can I reduce the additional space shown by blue arrows?

library(raster)

png(
  "Image.png",
  width     = 12,
  height    = 10,
  units     = "in",
  res       = 1200,
  point     = 10
)
 par(mar=c(1.5, 1.5, 1.5, 1.5), mfrow=c(3,4), oma = c(4, 4, 2, 2))

for (i in 1:12){

  month_col = i+2 #Since months start from 3 column
  r = rasterFromXYZ(testing[,c(2,1,month_col)])
  plot(r)

}

dev.off()


texts: Statistical bias arising out of only comparing US census block groups that have not changed (been split) over time?

texts:
Census Block Group ("CBG") study areas are notorious for changing over time. Even from 2009 to 2013, about 10-25% of the CBG's in South Carolina have not remained consistent over that time.

Seeing that CBG's provide the most geographic "resolution" for comparing data, I vastly prefer them over any other group. Trying to measure trends over time using them, however, can be difficult.

One way I'm trying to do this is by finding the CBG's that have not changed between 2009 and 2013 and measuring changes in just these areas. This way, even though I won't be able to get a complete picture of demographic changes over time, I will still be able to get a supposedly decent sample of the data changing over time within these tightly defined areas.

My concern is that maybe by only looking in these areas, I'm only looking at a certain type of area, namely, one that isn't changing enough to justify a "split" into multiple, smaller areas. Should I be worried about the statistical bias of only measuring the changes of CBG-defined areas that are not being redefined over time? If so, how much bias should I expect?


texts: Different in distance measurements in MapBasic/MapInfo

texts:
I am trying to make a custom ToolButton but I am having some problems getting correct measurements.

I have drawn a line (using the normal line tool) and set its length to 1000m. (Double-clicking it and measuring it with the ruler tool both show it to be exactly 1000.00m)

However when I use my tool and drag from one end of the line to the other it gives a distance of 999.8892m. I'm guessing it's some kind of projection issue but I'm using the same settings throughout the system... Where am I going wrong?

Include "mapbasic.def"
Include "icons.def"

Declare Sub Main
Declare Sub LineTool
Declare Sub ExitProg

Sub Main
  Alter ButtonPad "Drawing"
    Add
      Separator
      ToolButton
        Icon MI_ICON_ARROW_19
        HelpMsg "Click on map to draw a line"
        Cursor MI_CURSOR_CROSSHAIR
        DrawMode DM_CUSTOM_LINE
        Calling LineTool
    Show

  Alter Menu "Tools" Add
    "Close Line Tool" Calling ExitProg
End Sub

Sub ExitProg
  End Program
End Sub

Sub LineTool
  Set Distance Units "m"
  Set Coordsys Earth

  Dim x, y, x2, y2 as float

  x = CommandInfo(CMD_INFO_X)
  y = CommandInfo(CMD_INFO_Y)
  x2 = CommandInfo(CMD_INFO_X2)
  y2 = CommandInfo(CMD_INFO_Y2)
  print x & ", " & y & ", " & x2 & ", " & y2

  If x <> x2 and y <> y2 Then
    print Distance(x, y, x2, y2, "cm")
  Else
    print "That wasn't a line."
  End if
End Sub


texts: Displaying NoData transparent in QGIS diagrams

texts:
I am trying to display pie diagrams in QGIS 3.4. Lets say filling the whole circle would be 100%. If my data does not reach 100% in sum, I want the missing % to be empty or transparent.

Currently QGIS always fills the whole circle. Its expecting that the sum of my data is always 100%.

Here is an example:

Value_A;Value_B;Value_C;NoDataValue;Total
0,1;0,25;0,4;0,25;1
0,25;0,1;0,4;0,25;1
0;0;0,25;0,75;1
0,1;0;0;0,9;1
0,4;0,4;0;0,2;1

NoData is always 1-Sum(Value_A;Value_B;Value_C) and total is always 1.

Displaying Value_A, Value_B and Value_C in QGIS now looks like this:

So I was trying to add NoData and displaying it transparent. But transparent colors are not possible in QGIS diagrams. So it looks like this:

Basically I want to display the white areas transparent.

Is there some way (not necessarily by using transparent colors or diagram function) to achieve this task using QGIS 3.4? Using diagram size as %-sum is not an option in my case.


texts: .CSVT File no longer converting string variable to float value

texts: I was testing joining a dataset to shapefiles in Qgis. The text file that I used in my earlier testing as my csvt file worked perfectly and I was able to create gradient maps using my joined dataset. However, I have now moved the exact same dataset and the exact same csvt file over to another map in a different project. I joined the dataset to this map in exactly the same way with the csvt file in the same folder, just as before, but I am no longer able to create gradients using the data from the dataset. Specifically when I try to select the column to use for the gradient it should (and did) list 22 different variable options which are floats, now it does not list any as options. I have not changed anything in the text file, nor the dataset, so I have been unable to figure out why this worked fine earlier and no longer does.


texts: Extracting area weighted mean precipitation data not using a chart in Google Earth Engine

texts:
I am looking to just extract the area weighted mean value of hourly precipitation using the NLDAS-2 dataset in Google Earth Engine. I can do this using a chart, but it is very limiting considering it only produces 5,000 data points. I know there is a way to just export to table, but am having a challenging time doing so as I am a new user.
Right now I can get it to run with this code, but only 10 days at a time and I am looking to download 2002-2018.
// load in PRISM data and select dates
var PRECIP = ee.ImageCollection("NASA/NLDAS/FORA0125_H002")
                 .filterDate('2002-10-01','2002-11-30')

print(ui.Chart.image.series({
  imageCollection:PRECIP.select('total_precipitation'),
  region: Elv,
  reducer: ee.Reducer.mean(),
  scale:13875
}).setOptions({title:'hourly precip'}));


texts: Matching cities to values in .tiff file using ArcGIS Desktop?

texts:
I am new to ArcGIS Desktop.

I have a list of cities and their longitude/latitude (as a .csv). I have another .tiff file which contains a map of the world records values on a 5' x 5' grid cell (this is the data: https://ozak.github.io/Caloric-Suitability-Index/).

I am now trying to match my cities to the average value of the .tiff data in a 50km radius around the city. I would like to export the data as a .csv or something similar with an added column for the average value.

I've pulled both datasets up in ArcMap but I cannot figure out how to overlay, calculate the average and then export the data again.

Could someone point me in the right direction?


texts:
When I export to PDF, I get some black lines either vertical or horizontal.
It's exported in 300 DPI and I have used QGIS 2.8.


texts: Model Builder not showing feature layer FID

texts:
In the middle of a routine I'm making on model builder I have a generate near table function beetween two feature of points. The results give me a table with only the coordinates of the points and their FIDs of every combination. From these coordinates I make lines using the XY to line and keep all the FIDs on this new shapefile.

What I want to do is use the join field function to recover all the information in the attribute table form both points of each line using the FID from both tables (line feature and point feature) as the common attribute. However the model builder doesn't recognize that there is a FID field on the points attribute table. If I do the whole routine myself without model builder it works, but on model builder I can't work around it.

Anyway I can fix it?


texts: Adding compression to a NetCDF file using xarray

texts:
I have a raster that I am trying to compress and convert to a NetCDF format with compression level = 9 using the xarray package. I assume that the compression is added using the encoding parameter as a dict, but I am not quite sure how I understand what to add completely here:
 f = directory + "/D_Passaic_F02_NBR_E0001_WGS84_comp"
 t = xarray.open_rasterio(f)
 encode = {'zlib': True, 'complevel': 9}
 t.to_netcdf(output_dir+"/Test2.nc", encoding=encode)

I basically get an error: KeyError: 'zlib', but I am not sure what I am suppose to use here. Suggestions?
The data array shows a single band, and an x, and y variable like this:
 <xarray.DataArray (band: 1, y: 9635, x: 14564)>

        or, in more detail:

 <bound method ImplementsArrayReduce._reduce_method.<locals>.wrapped_func of
 <xarray.DataArray (band: 1, y: 9635, x: 14564)>
 [140324140 values with dtype=float32]
 Coordinates:
   * band     (band) int32 1
   * y        (y) float64 41.06 41.06 41.06 41.06 ... 40.74 40.74 40.74 40.74
   * x        (x) float64 -74.45 -74.45 -74.45 -74.45 ... -73.97 -73.97 -73.97
 Attributes:
     transform:                 (3.2670488250568696e-05, 0.0, -74.447024371179...
     crs:                       +init=epsg:4326
     res:                       (3.2670488250568696e-05, 3.2670488250568696e-05)
     is_tiled:                  1
     nodatavals:                (-9999.0,)
     scales:                    (1.0,)
     offsets:                   (0.0,)
     AREA_OR_POINT:             Area
     HISTOGRAM:                 9090|9307|9097|9209|8729|8864|8744|8864|9181|9...
     TIFFTAG_ARTIST:            HEC-RAS
     TIFFTAG_IMAGEDESCRIPTION:  Depth (Max)>


texts: Editing scripts outside (online) Code Editor?

texts:
While reading EE docs https://developers.google.com/earth-engine/playground, something called my attention:

"Repositories can be accessed using Git, so you can manage and edit your scripts outside the Code Editor..."

I use EE online editor and I'm aware of possibilities git bring (I think). But I'm wondering in which circumstances it might be useful to be able to "...edit your scripts outside the Code Editor..."? I guess Code Editor mean the online one.

I understand, outside the online code editor the scripts might be edited as any file. But I cannot run/test them, at least not before uploading them back to EE. So I'll be doing some kind of blind changes. Perhaps I'm missing something.


texts: I'm confused as to the range of brightness values associated with Landsat 8 Bands. I know that for Landsat 7 and below it has always been 0 - 255 but some of the literature indicates that Landsat 8 values range from  0 - 4096 since it's 12 - Bit. Is this a correct?


texts: Using MapServer with ESRI enterprise geodatabase on Oracle

texts:
I've been looking into MapServer, but can't find anything indicating that it can connect to an ESRI enterprise geodatabase (version 10.3 or 10.4) for vector data. There is ArcSDE support but it's limited to version 9.2.

Is there any known path to serving vector data from an enterprise geodatabase on Oracle that uses the ESRI ST_GEOMETRY format?


texts: Drop multiple fields in graphic modeller qgis 3.10.14 not working

texts:
I am trying to use the graphical modeler to delete multiple fields in QGIS 3.10.14.

When i try to use the Drop Fields algorithm to delete multiple fields from the required layer, it deletes only the first field and not all the specified fields.

From the console result it states that the other fields not deleted don't exist.

I tried using different delimiters to separate each value i.e. comma, semicolon. It did not work either.
However if I run the Drop Fields multiple times, it removes all field in each of the operations.

So will it be possible to delete the fields in one operation?


texts: Converting *.tif file into headerless (flat) binary file?

texts:
I am using NDVI files (value range -1 to +1) created from Landsat CDR images to produce a timeseries in TIMESAT software.

My  .tif file has following characteristics.
File type  TIFF,Data Type  Double Precision, Byte Order Host (Intel), Interleave BSQ

As per the requirement of TIMESAT, the image should be headerless(flat) binary file, and values can be stored as 8 bit integer(0-255), 16 bit signed integer(-32767 to +32768) or 32 bit real.

Now i need to change my data type into any of these 3 format and .tif file into headerless binary files.

Any idea?


texts: Identifying road intersections using PostGIS

texts:
I'm trying to identify where roads intersect each other, and to make a point at this intersection, with the number of roads that form the intersection listed.

I was wondering if there was some way of using ST_NumPoints to achieve this but I can't quite figure out what I should be doing. What I have done is to create a table of points where the lines intersect using the following code:

CREATE TABLE test_points as
SELECT
    ST_Intersection(a.geom, b.geom),
    a.gid
FROM
    roads as a,
    roads as b
WHERE
    ST_Touches(a.geom, b.geom);

If I run this on a sample of roads I get the following grid of points (the roads are shown for illustration):

If I inspect one of the points, I see that there are many points stacked on top of each other:

The GID here is the road ID, but I don't understand why there are some many points. I can understand 4 points being counted for a central road intersection, but there are 12 points listed here. Is there a better way to perform this calculation in PostGIS?


texts: Anyway to integrate QGIS in Visual Studio Enterprise/Community?

texts:
I would like to ask if it possible to integrate QGIS in Visual Studio (either enterprise or Community)?
You see I am developing for  a python based application which request JSON data from an API. The Spatial analyses and other stuff happens in QGIS more specifically the PYQT part of the GIS API.
Is this possible? Currently I am trying to work with three different Python Envoriments: Iron Python 2.7 for a GUI, Python 3.8 for conversion of the JSON file and PyQT for conversion of the file in spatial terms.
Yes I am aware of QT designer but considering I am more familiar with Visual Studio and with Iron Python for designing a From, it is the best option for me.
Keep Note Iron Python is just being used to develop a WPF/Form to return a reponse url, I have heard some storys that PYQT does not work with PYQT. So I am uncertain of certain things.
So possible or not?


texts: Has QGis 3.4 some limitation for filenumber to load shp files?

texts: I want to load over 600 shp files at once in qgis 3.4. Is there a limitation on the number of shapefiles that can be added?


texts: ST_Buffer with segment AND endcap parameters

texts:
Making a QGis plugin that buffers a line and I want to declare a high number of segments and a flat endcap.

ST_Buffer(geom, 40, 25, 'endcap=flat')

This doesn't work with both parameters entered, it seems I can only choose one or the other.

How can I, if I can, choose both?

Current versions of both pieces of software.


texts: "Too many points" when using gdal2tiles.py in EPSG:4326

texts:
I want to put a custom map onto cartodb, which, I believes, requires EPSG 4326. When I try making tiles with gdal2tiles, using a 4326 georeferenced image, I get the following error:

ERROR 1: Too many points (10201 out of 10201) failed to transform,
unable to compute output bounds.

If I try to use other projections, gdal2tiles will work fine but they show up in the completely wrong place on cartodb.


texts: Automate daily reclassification in ArcGIS

texts:
I have a sample raster of my area which is unclassified, and I have a MS-Excel file with the first column for the day of the year and the second column showing the corresponding daily value in %. This % value will vary from day to day, hence for 365 days, I will have 365 rows for each day & its respective % value. Normally, I would use the Reclassify tool in ArcGIS to manually define the breaks (by referring to the MS-Excel %value column) for instance-Please note that the old values column in the screen-shot were taken from the MS-Excel file % value column for that particular day:

Thus, for each day in a year, I will have to do this manually in ArcGIS, so is there a way to automate this in R (maybe the raster package!) or in Python? Please note that the Excel file is not gridded!
UPDATE
The % values of each day will vary from each other, so in the script (R or Python), I will have to manually enter the ranges using a loop and If-else for each day. Basically, New values 1 i.e.Old values 0-39 remain constant for each day, so this break value of 39 will stay the same for each day. It is the New values 2 that change in accordance with the % value in the Excel file (96% for the day seen in the screenshot) of that day. And, New value 3 is the remainder %. But, If I could read the MS-Excel file into the script that basically would automate the break values for New values 2. I hope this further clarifies the question!


texts: OpenLayers WFS protocol query after vector layer creation

texts:
I want to be able to create a layer on init, then load features from a WFS protocol afterward. I'm having trouble doing this.

Using WFS protocol in vector constructor works fine and is not my problem:
http://dev.openlayers.org/releases/OpenLayers-2.10/examples/wfs-filter.js

What i want to do is Add a layer like such:

var layer = new OpenLayers.Layer.Vector("Results", {
    projection: new OpenLayers.Projection("EPSG:4326")
});

Then Query it using the WFS Protocol.... this is what i can't get to work... what i have so far is:

var protocol = new OpenLayers.Protocol.WFS({
    version: "1.1.0",
    url: MapBase.Configuration.GeoServer.WFS_Service.EndPointURL,
    featurePrefix: MapBase.Configuration.GeoServer.WFS_Service.featurePrefix,
    featureType: MapBase.Configuration.GeoServer.WFS_Service.featureType,
    featureNS: MapBase.Configuration.GeoServer.WFS_Service.featureNS,
    geometryName: MapBase.Configuration.GeoServer.WFS_Service.geometryName
});

var response = protocol.read({
    maxFeatures: Search.Properties.maxResults,
    callback: function (resp) {
        console.log(resp);
        layer.addFeatures(resp.features);
    }
});

resp.features gets populated fine. My main problem is I can't find where to specify my filter and  projection as i can in the Vector constructor. Even though i specified the projection in the vector constructor, when i addFeatures from my WFS response, they come back in a different projection because they end up really small spooning the west coast of Africa. Same thing with the filter. I don't see that WFS protocol alone supports applying a filter? It only seems to work in conjunction with a vector layer constructor.

This seems wrong though because the WFS protocol as described by the OGC support the filter and when openlayers submits it's WFS request via the vector constructor, it does indeed include it just fine. There should be a way to specify projection and filter in the WFS protocol alone?

EDIT #1:

What i have now:

var protocol = new OpenLayers.Protocol.WFS({
    version: "1.1.0",
    url: MapBase.Configuration.GeoServer.WFS_Service.EndPointURL,
    featurePrefix: MapBase.Configuration.GeoServer.WFS_Service.featurePrefix,
    featureType: MapBase.Configuration.GeoServer.WFS_Service.featureType,
    featureNS: MapBase.Configuration.GeoServer.WFS_Service.featureNS,
    geometryName: MapBase.Configuration.GeoServer.WFS_Service.geometryName,
    srsName: "EPSG:4326"
//    ,
//    defaultFilter: myFilter
});

protocol.mergeWithDefaultFilter(Search.GetWFSSpatialFilters());


texts: Removing Google maps base layer from Google Maps API

texts: I need to remove the base layer of google maps dynamically. User can add or remove the base layer.Because I am using more layer on top of it.Is there any function in google maps api to remove the basemap and leave it as a blank.


texts: Conversion of Fbx or Obj to Cesium Tiling

texts:
I have a an OBJ File and a FBX File of very large size (OBJ = 25GB & FBX = 6GB)

I'm looking for a way to convert either of them into Cesium Tile Format. I have tried using Pix4D software which is capable of converting it but it reduces the model quality significantly. So kindly suggest something that doesn't affect the quality much.


texts:
Is it possible to get categorization of a gpx file in Google my maps by elevation. Something like this:

When I import gpx file I am able to get displayed elevation attribute of the whole route and from first and last point. I guess its because of automatic calculation of the software with which I recorded map. I guess that this photo above, which was done with online software http://maplorer.com/view_gpx.html probably does some interpolation.

I am explicitly looking for a manual solutions, not coded one.


texts: Splitting field into several fields based on strings separated by space in QGIS?

texts:
I need to create dichotomic fields based on a single field which shows multiple choice answers in the form of strings separated by a space. What I want to do is this:

The number of choices that will produce the new fields is known, but I have 27 fields with data from multiple choice questions in several layers, I know how to do this manually, but it would take a very long time.

update:
Here is a screenshot of one of my field of interest, there are blank spaces, there are no "0" because the multiple-choice options begin with "1". When in a group this field is not filled (because it depends on previous selections) then it becomes blank.


texts: Geometry type issue with st_intersection() between polylines and polygons

texts:
I clip a polylines vector layer with a polygon vector layer, containing one single polygon. I need my polylines to be clipped at the polygon's boundaries.

CREATE TABLE test AS (SELECT a.gid, st_intersection(a.geom, b.geom) FROM lines a, polygon b WHERE st_intersects(a.geom, b.geom));

I can't visualize the test table in QGIS because the type of geom column seems not to be recognized (see image below, which is DB Manager GUI). Thus I tried to specify it within pgAdmin:

ALTER TABLE test ALTER COLUMN geom TYPE geometry(MultiLineString) USING geom::geometry(MultiLineString);

This issues the following error message:

  ERROR:  Geometry type (LineString) does not match column type (MultiLineString)

--

EDIT :

I tried to keep the geometry column from the polylines layer together with the st_intersection() function.

CREATE TABLE test AS (SELECT a.gid, a.geom, st_intersection(a.geom, b.geom) FROM lines a, polygon b WHERE st_intersects(a.geom, b.geom));

This allows to display the table within QGIS but as there are 2 geometry columns produced, the clipping is not displayed and QGIS considers that 2 tables were created.

I need to make a clean clipping of my polylines layer resulting in a new table with a geometry columns having geometry(MultiLineString) type and displaying well in DB Manager.


texts: Is there an automatic backup script for ArcPAD 10 when you are out in the field?


texts: pyinstaller package of pyqgis 3 standalone program failed import qgis.core

texts:
I have a working standalone pyqgis program with mapping ability, it was in python 2.7 and qgis 2.x library. I used pyinstaller to package it up. Now I have migrated this application to python 3.6 and qgis 3.2.1-Bonn and pyqt 5.9. It is working fine in the PyCharm environment, the map is working, basically, qgis imports are working, such as:

from qgis.core import QgsApplication, QgsVectorLayer

I tried to package it up with pyinstaller 3.3.1 (as I have done successfully for the python 2.7 version of the program), it can get all of the qgis pyd/dlls into the package and the Qt GUI is running, it just keeps failing the qgis import statements such as above. The whole package is gathered into the C:\dist\qgis_test\ folder that contains the following qgis libs:

qgis._core.pyd
qgis._gui.pyd
qgis_core.dll
qgis_gui.dll
qgis_native.dll

When running the packaged executable file, it always bombs at the point of importing modules from the qgis.core or qgis.gui packages as shown below:

PYIMOD03: fullname=qgis._gui
PYIMOD03: filename=C:\dist\qgis_test\qgis._gui.pyd
PYIMOD03: fullname=qgis._core
PYIMOD03: filename=C:\dist\qgis_test\qgis._core.pyd
Traceback (most recent call last):
  File "qgis_test.py", line 4, in <module>
    from qgis.core import QgsApplication, QgsVectorLayer, QgsProject
  File "C:\OSGeo4W64\apps\Python36\lib\site-packages\PyInstaller\loader\pyimod03_importers.py", line 631, in exec_module
    exec(bytecode, module.__dict__)
  File "qgis\core\__init__.py", line 27, in <module>
  File "C:\OSGeo4W64\apps\Python36\lib\site-packages\PyInstaller\loader\pyimod03_importers.py", line 716, in load_module
    module = loader.load_module(fullname)
ImportError: DLL load failed: The specified procedure could not be found.

The thing is that all of the PyQt5 libraries are all imported without error.

I've already set the following environmental variables to the package root directory (because that's where qgis_core.dll and qgis_gui.dll are located) like I did for the python 2 package:

QGIS, QGIS_PREFIX_PATH, PATH, PYTHONPATH, QT_PLUGIN_PATH

I pretty much searched and read up all python3 qgis related threads here and via google, but not much information is out there. Also, I noticed that the Qt plugins are under PyQt5\Qt\plugins within the package directory, I also set the environment variables accordingly to that, but still no luck. It keeps failing on those qgis import statements.


texts: Calculating Total Line Length of a Network with Shapely: AssertionError

texts:
I would like to calculate the total length of a line network, similar to what is done in this QGIS tutorial.

Here is what I have:

import fiona
from shapely.geometry import LineString

length = 0.0
with fiona.open('/Users/felix/data/centerline.shp', 'r') as shp:
    print(shp.crs)
    for rec in shp:
        line = LineString(rec['geometry']['coordinates'])
        length += line.length

While this method works fine for some shapefiles, for others it produces this error:

{'datum': 'NAD27', 'proj': 'lcc', 'y_0': 0, 'lat_2': 35.46666666666667, 'lat_0': 33.5, 'no_defs': True, 'x_0': 609601.2192024384, 'lon_0': -118, 'lat_1': 34.03333333333333, 'units': 'us-ft'}
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 200, in geos_linestring_from_py
    array = ob.__array_interface__
AttributeError: 'list' object has no attribute '__array_interface__'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/felix/Documents/Scripts/python/county_workflow_automation/scratch.py", line 11, in <module>
    line = LineString(rec['geometry']['coordinates'])
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 49, in __init__
    self._set_coords(coordinates)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 94, in _set_coords
    self._geom, self._ndim = geos_linestring_from_py(coordinates)
  File "/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/shapely/geometry/linestring.py", line 271, in geos_linestring_from_py
    assert n == 2 or n == 3
AssertionError

Can somebody help me interpret this error? Or, if possible, is there a more direct way to accomplish this task?


texts: Advanced search terms in the Semi-automatic Classification plugin

texts:
Background:
At QGIS 3.18.1, I can successfully download Sentinel Level-1C imagery by entering the granule that covers my area of interest (10UGU) in the Advanced search box.
Now I would like to download Level-2A imagery for the same granule.
The Semi-Automatic Classification Plugin documentation states on page 48:

"... some services (e.g. Sentinel-2) allows for entering search
parameters; for instance, it is possible to enter the name of a
granule (e.g. 33TTG) without defining a search area in order to get
the results of that granule; Sentinel-2 level 2A can be filtered
entering S2A_MSIL2A;*" [my emphasis].

What I want to do now is have two search parameters in the Advanced search box, one for the granule (10UGU), and the other for the Level-2A imagery (S2A_MSIL2A*;).
Problem:
Whenever I add the two parameters in the Advanced search box, and then click the Find button to search for available imagery, SCP responds with the statement Information [23] Define a search area as shown here:

I have tried adding the two parameters in every combination that I can think of: placing one before the other, separating them with commas, semicolons, spaces, etc.  In every case, SCP responds with the error shown in the screenshot above.
Question:
How do I add multiple search parameters to the Advanced search box?


texts: How to batch geoprocess using iteration with tool parameter values in table

texts:
As part of my work at a non-profit, I am preparing a map for an upcoming legislative session. GIS is only a small portion of my job. I am comfortable with ArcGIS but have no experience with Python (I do know Java, Visual Basic, and HTML and have a basic understanding of object oriented coding). Here is the issue (that I think can be relatively easily done in Python).

I want to generate a certain number of random points in roughly 1,200 ZCTAs. The number of points corresponds to the client records for that ZCTA. I have a table with the ID for of each ZCTA in the state and the number of client records for that ZCTA. I also have a shp file of all ZCTAs in the state. Given the number of ZCTAs, doing this manually will take too much time.

The table I have has N rows and 2 columns, ID and PointNum. The data table of the ZCTA shp data has one field, GeoID. The table can be joined to the shp with the ID field.

The idea is to generate the correct number of random points in each ZCTA by iterating through the table and setting the parameter of the random point tool to the value of PointNum for the given ZCTA. Here is a conceptual model for how I could do this.

WHILE x < N
Select feature in ZCTA shp using table ID field
Random point tool with extent same as ZCTA shp and number of points = PointNum
Clear Selection
N + 1

I have no idea how to implement this (or even if my proposed way makes sense). I seek an example of a python script (with some annotation).


texts: Does QGIS automatically use spatialite spatial indexes?

texts: When you write a query in spatialite you have to manually join your query to the spatial index for that index to be used to speed up the query.  When QGIS displays data does it take advantage of any spatial indexes as you pan and zoom around your data?


texts:
I'm trying to add a vector tile layer in my app using OpenLayers 4, but I'm getting the error as follows:
ol.js:753 Uncaught TypeError: a.replace is not a function

Here is my code to add the vector tile layer,
var vectorSource = new ol.layer.VectorTile({
        declutter: true,
        source: new ol.source.VectorTile({
            format: new ol.format.GeoJSON(),
            url: function (extent) {
                return 'http://localhost:8080/geoserver/myworkspace/wms?service=WMS&version=1.1.0&request=GetMap&layers=myworkspace:bdms_building&styles=&bbox=227227.0625,2664047.75,300691.15625,2745836.25&width=689&height=768&srs=EPSG:32640&format=application/x-protobuf;type=mapbox-vector&' +
                    'bbox=' + extent.join(',') + ',EPSG:32640';
            }
        }),

    });

        proj4.defs("EPSG:32640", "+proj=utm +zone=40 +datum=WGS84 +units=m +no_defs");
        var map = new ol.Map({
        target: 'map',
        controls: ol.control.defaults({
            attributionOptions: {
                collapsible: false
            }
        }).extend([
            new ol.control.ZoomToExtent({
                extent: [
                    227227.0625, 2664047.75,
                    300691.15625, 2745836.25
                ]
            })
        ]),
        interactions: [
            interactionSelectPointerMove,
            new ol.interaction.MouseWheelZoom(),
            new ol.interaction.DragPan()
        ],
        layers: [
            vectorSource
        ],
        view: new ol.View({
            center: [227227.0625, 2664047.75],
            projection: 'EPSG:32640',
            maxZoom: 19,
            zoom: 13
        })
    });


texts: Create catchments/subbasins for specific nodes

texts:
I have an XPSWMM model of a neighborhood storm drainage system. Nearly all of the drainage is above ground through a network of ditches and culverts.

I still need to create catchments for the upstream points for each culvert so I can route runoff through the network. There are around 150 culverts, and the topography is quite complex, and doing this manually is quite time-consuming and complicated.

Is it possible to create somewhat of a network of subcatchments to snap to these points in QGIS?


texts: Using a Python Script Input parameter as an output

texts:
I have a Python script that takes several input parameters, one of which is the name of an output folder which will hold processed data (outFolder). I want to be able to access that outFolder variable displayed as an output (outFolder--> Process-->outFolder) as well once the process has finished and "outputted" the results into the outFolder. Basically, I want to use that same filled outFolder as an input for another process in ModelBuilder. How do I do this?

#Get the input feature class, optional fields and the output filename
input_FC = sys.argv[1]
inField = sys.argv[2]
theFName = sys.argv[3]
outFolder = sys.argv[4]


texts: Duplicating raster layer using QGIS Modeler

texts:
I have a repetitive process where I duplicate a raster layer twice so I have 3 instances and I can display hillshade, elevation, and contour styles - I'm trying to create a model to run this process for me, but can't seem to find a duplicate raster algorithm.
Am I missing something or is this not possible with model designer?


texts: Showing only legend items that are actually displayed within current map extent using ArcGIS API for JavaScript?

texts:
I'm using TOC/Legend widget as a default legend panel in my ArcGIS based JavaScript API web app.

It works great, beside now i want to filter this legend panel based on users current map extent. The reason is, currently i have some 30 to 40 items in this TOC pane and it will better to filter them out based on the location user is focusing on. So, if any extension or widget someone already using or modifying the existing TOC will help.

There is already an ArcScript present for ArcGIS Desktop called  Legend Limiter, but that is not for JavaScript API. Also a similar question is asked as Showing only legend items that are actually displayed on map in ArcGIS Desktop? but again, that doesn't address the JavaScript API issue.


texts: How not to Copy Features when too many would be output?

texts:
This question has already been posed to the ArcGIS Geoprocessing forum but I have had no responses there.

I have a 3.5 million polygon layer (cadastre) stored in Oracle via ArcSDE. I use a polygon feature chosen by an end user from a second layer to copy out the cadastral polygons that intersect it (via Copy Features). I can get this to work fine for most polygons, and get the few seconds to few minutes response I need, by using the polygon the user chooses to set the geoprocessing extent (arcpy.env.extent) prior to executing Copy Features.

However, for a few polygon features that the user might choose (think petroleum or gas pipeline corridor), the number of polygons that would be copied out can number up to about 200,000. These can be copied in about 15 mins, but rather than get all of these back so that I can write them into a very bloated PDF report, I would like to save time and pages by having the option to set a maximum number of features that can be copied out i.e. be able to set a MaxRecordCount = 500 that would stop the Copy Features once 500 had been copied.

Is there a way to tell Copy Features to only copy up to a maximum number of features?

Alternatively, is there a fast way to get a count of how many features in an ArcSDE layer are within the current geoprocessing extent?

UPDATE: No real progress made but the results and code below may illustrate better how long it takes to use SelectLayerByLocation and GetCount (at least at ArcGIS 10.0) against the 3.5 million polygon ArcSDE layer even when the tenement chosen only overlaps 9,000 polygons (original example takes 10x as long).

This is the output:

Copying one tenement to in_memory takes 2.87037144971 seconds
Setting arcpy.env.extent to that one tenement (in_memory)
SelectLayerByLocation takes 64.4129756363 seconds to intersect cadastre with this one tenement (in_memory)
8961 cadastral parcels found in Search Area
GetCount takes 206.410587153 seconds to come up with that number

This is the code:

import arcpy
import time

mxd = arcpy.mapping.MapDocument("C:\\MXDS\\CustRep\\Maps\\Search_layers_V1.mxd")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
searchAreaLayer = arcpy.mapping.ListLayers(mxd,"Exploration Permits Petroleum (EPP)",df)[0]
searchAreaLayer.definitionQuery = "TENURENO = 644"
start = time.clock()
arcpy.CopyFeatures_management(searchAreaLayer,"in_memory/SearchArea")
elapsed = (time.clock() - start)
print "Copying one tenement to in_memory takes " + str(elapsed) + " seconds"

print "Setting arcpy.env.extent to that one tenement (in_memory)"
arcpy.env.extent = "in_memory/SearchArea"

cadastreLayer = arcpy.mapping.ListLayers(mxd,"Lot and Plan",df)[0]

arcpy.MakeFeatureLayer_management("in_memory/SearchArea","SearchArea")

start = time.clock()
arcpy.SelectLayerByLocation_management(cadastreLayer,"INTERSECT","SearchArea")
elapsed = (time.clock() - start)
print "SelectLayerByLocation takes " + str(elapsed) + " seconds to intersect cadastre with this one tenement (in_memory)"

start = time.clock()
result = arcpy.GetCount_management(cadastreLayer)
numRowsInSearchArea = int(result.getOutput(0))
print str(numRowsInSearchArea) + " cadastral parcels found in Search Area"
elapsed = (time.clock() - start)
print "GetCount takes " + str(elapsed) + " seconds to come up with that number"


texts: How do I know the IDs of the multiple links that are involved in the shortest distance between any two points in a Network?

texts:
I don't just want to find the shortest distance between any points in the Network Dataset but find which links (the IDs) are involved in the shortest path.

At the end I want a table that includes the origin, destination, and the links involved in the shortest distance between origin and destination. I have around 80,000 rows in the Origin-Destination Table.


texts:
Trying to generate a vector-based PDF map, which includes transparencies, for quality purposes.  However, I'm finding that transparencies get rasterized instead and is lower quality.

Is there any solution to this other than exporting the PDF to Illustrator and manually applying transparencies there?


texts: Satellites with higher revisit frequency

texts:
Sun-Synchronous orbit has a 16 day revisit. Geo synchronous though further away, observe the same spot on earth, and can provide data with higher frequency.
Most earth observation satellites that share data publically I know are in sun-synchronous orbit except GOES (east and west), which provides data every 5 mins.

I am looking for satellites that can provide data with higher frequency. Are there other satellites which provide data at higher frequency? (We don't care about what location they are observing on earth)


texts: Download Solar Insolation Data from NASA in R

texts:
I need to get solar insolation data from different parts of the world. I have several latitude and longitude coordinates and I would like to get a mean solar insolation value across different years.

I found this data from NASA:

https://neo.sci.gsfc.nasa.gov/view.php?datasetId=CERES_INSOL_M

From that website you can download the GEOTiff files, but I couldn't find a way to download several files (say monthly) across different years to get a mean estimate.

Is there a way to do it? Or are there any other resources that I could use to get this data (preferably in R)?


texts: Creating WHERE clause from list using ArcPy

texts:
I am creating a new shapefile (a subset of a current shapefile) such that it has all the features whose values in a field ("fclass") is not contained in the field_names list.

field_names = ['bridleway','cycleway','footway','living_street','pedestrian','residential','tertiary_link','service','unclassified','steps','path','tertiary']

def get_where_clause(field_values):
   where_clause = ""
   for value in field_values:
      where_clause = where_clause + """"fclass" <> '{0}'""".format(value) + " AND "
   return where_clause[:-5] #removing the last "AND"

arcpy.MakeFeatureLayer_management(path, feature_layer, get_where_clause(field_names)

Problem:
The output somehow selects all the features and exports everything.


texts: Set schema in ogr2ogr from GML to pgdump

texts:
My ogr2ogr tranformation from a GML to a pgdump sql file always by default uses the schema public. I would like to use a different schema. It seems to be possible: How to assign a schema to ogr2ogr?. However in my dump file every line still states the public schema.

This is the code that does run but has the output with still the public schema:

ogr2ogr -f PGDump file.sql gi*.xml -oo GML_ATTRIBUTES_TO_OGR_FIELDS=YES -lco SCHEMA=klicbestanden_in


texts:
I'm new to the open-source GIS world, having abandoned ArcGIS due to costs, and am looking for an OS analog to the Collector app. What I'm seeking to do is collect field data with an Android cell phone (namely points, tracks, polygons to which I can assign some attributes) and import it into QGIS, like I'd have done with Collector and Arc. I've toyed with NextGIS, GeoODK Collect, and QField but they seem crazy confusing to a noob like myself.

Any hints or recommended starting points?

I can get basic GPX data from any old mapping app but I'm looking for something more robust and scalable.

UPDATE: I've gotten QField to the point where I can create an informative map in QGIS and upload it to my phone; the issue now is figuring out how to record GPS tracks and the like. I'm finding the clunky interface (and online documentation) a bit lacking. I gather I can simply draw polygons and lines, but that's not nearly as helpful as being able to log my GPS track. I need this functionality to map logging roads/trails. Any ideas?


texts:
I have a csv file containing information for loads of lines (dates, road names, condition information, etc etc) but which have no geographical information (other than the association of road name)

Is it possible to load the file in QGIS, and draw the polyline which refers to each row of data, rather than have to re-enter the data for each polyline drawn?


texts: Extracting boundary of raster into polygon

texts:
I have elevation data downloaded from nationalmap.gov in ArcGrid format (1/3 arc-second DEM). My main goal is to obtain the boundary / extent of each raster range that I have created by simply extracting by attributes. For example, I have a raster within the range of 3 meters and 10 meters (see attached image) and I would like this to be a shapefile and preferably a much smaller file size. It is currently 15gb uncompressed, and I cannot simply convert to shapefile as it is not given as an option when I try to use that tool. I have also tried Raster Domain via 3D analyst, but the result is still far too large and also unable to be converted to a polygon.

What would be the best approach for me to take to reach my objectives?

I have a lot of these ranges that I need to convert and then utilize in a leaflet interface.


texts: Texts are smaller when I open project file by another computer

texts: I have made a project with text on map layout. But when I opened that project by another computer, I found the texts are smaller (minus 1 size). How do I solve this problem?


texts: Implementing new Snapping in ArcGIS Engine C# Application?

texts:
As I can find, there are two kinds of snapping in ArcGIS 10.2.2:

1- Classic snapping.

2- New ArcGIS Snapping, which doesn't have any specific name as I know.

I can't implement the New ArcGIS snapping in an ArcGIS Engine application. As I have read, it is completely available in ArcGIS Engine library but I can't make it work. The following is the code I'm working on :

IHookHelper2 m_hookHelper2 = new HookHelperClass
 {
       Hook = mapMain.Object
 } as IHookHelper2;
 IExtensionManager extensionManager = m_hookHelper2.ExtensionManager;

 if (extensionManager != null)
 {
      UID guid = new UIDClass();
      guid.Value = "{E07B4C52-C894-4558-B8D4-D4050018D1DA}";
      IExtension extension = extensionManager.FindExtension(guid);
     ISnappingEnvironment m_SnappingEnvironment = extension as
                ISnappingEnvironment;
     m_SnappingEnvironment.SnappingType =
                esriSnappingType.esriSnappingTypeEdge |
                esriSnappingType.esriSnappingTypeEndpoint;
     m_SnappingEnvironment.Enabled = true;
     m_SnappingEnvironment.ShowSnapTips = true;
     m_SnappingEnvironment.SnapTipType =
                esriSnappingTipType.esriSnappingTipLayerName;
     ISnappingFeedback m_SnappingFeedback = new SnappingFeedbackClass();
     m_SnappingFeedback.Initialize(null, m_SnappingEnvironment, true);

     if (m_SnappingFeedback != null)
     {
           m_SnappingFeedback.Refresh(0);
     }
 }

I should mention again that I want to implement New ArcGIS Snapping not the classic one.


texts: How to select all points of a polygon in Spatialite

texts:
I need to select all vertex coordinates from a polygon table in Spatialite DB. In PostGIS, I'd normally do something like

SELECT
  st_x((st_dumppoints(poly_geom)).geom),
  st_y((st_dumppoints(poly_geom)).geom)
FROM
  my_poly_table;

But there is no st_dumppoints function in Spatialite.
How can I get a table with every vertex of every feature?


texts:
I'm trying to snap two Linestrings shapefiles using Shapely/Geopandas. I tried :

import geopandas as gpd
from shapely.geometry import *
from shapely.ops import snap

lines1 = gpd.GeoDataFrame.from_file('lines1.shp')
lines1.to_crs({'init': 'epsg:2227'})
lines2 =gpd.GeoDataFrame.from_file('lines2.shp')
lines2.to_crs({'init': 'epsg:2227'})
res = snap(lines1, lines2, 10)
res.to_file('result.shp', driver="ESRI Shapefile")

But I got this error :

Traceback (most recent call last):
File "C:/Users/Documents/script.py", line 14, in <module>
    res = snap(lines1, lines2, 10)
File "C:\Program Files (x86)\Python27\lib\site-packages\shapely\ops.py", line 303, in snap
    return(geom_factory(lgeos.methods['snap'](g1._geom, g2._geom, tolerance)))
File "C:\Program Files (x86)\Python27\lib\site-packages\pandas\core\generic.py", line 2672, in __getattr__
    return object.__getattribute__(self, name)
AttributeError: 'GeoDataFrame' object has no attribute '_geom'

Process finished with exit code 1.

Any idea on what the problem is?

EDIT :

Updated code based on @caiohamamura's answer :

import geopandas as gpd
from shapely.geometry import *
from shapely.ops import snap

lines1 = gpd.GeoDataFrame.from_file('lines1.shp')
lines1 = lines1.to_crs({'init': 'epsg:2227'})
lines2 = gpd.GeoDataFrame.from_file('lines2.shp')
lines2 = lines2.to_crs({'init': 'epsg:2227'})
res = lines1
lines2_union = lines2.geometry.unary_union
res.geometry = res.geometry.apply(lambda x: snap(x, lines2_union, 10))
res.to_file('result.shp', driver="ESRI Shapefile")

Now I have this error :

Process finished with exit code -1073741819 (0xC0000005)

EDIT 2 :

I no longer get the above error. The script runs well, although the snapping doesn't seem to work that well. Here's an example of the result :

The blue line is the result of the snapping

The black one is lines2

The red one is lines1

I tried changing the tolerance but it seems that if I give a value that's > to 15 feets I might get a wrong output.
The lines do not get snapped totally to each other. In some places they do, in others they don't as demonstrated on the screenshot above.


texts:
Preface: I'm a total GIS noob. Flame away. I'm learning GIS in my free-time and with free resources.  I realize that is likely going to force me on the learn-fast-or-giveup track, but here it goes...

My first project: Create a laser cut 9-layer topographic map of the 7.5 minute area where I work.

My Plan:

1) Download SDTS DLG Hypsography data from USGS.

2) Unzip using Peazip.

3) Import into GRASS GIS.

4) Remove/edit unnecessary topo lines (I was able to do it with dlgv32 available at USGS website, but I can't export and don't want to pay $500 for a full version of global mapper).

5) Export to DXF or DWG

6) Open in Corel and create layers

7) Laser cut layers, glue them together.

I've gotten through steps 1, 2, and 3 with minimal frustration, although peazip throws an error on extract that says "data after archive", which I hope has no effect on the usability of my data. I'm at a loss on step 4.

That data set I'm working with is DLG STDS data in large format (7.5 minute quadrangles).  Specifically, the data is from: http://dds.cr.usgs.gov/pub/data/DLG/LARGE_SCALE/I/indian_hills_CO/hypsography/1/

I can get the DDF files to import into GRASS and I can view the topomap in the display; however I'm having trouble manipulating what I'd like to see. Basically I want to display only the major topolines (for example, 6000', 6400', etc.).  When I open this data set in Global Mapper, this is very easy to do.  I imagine I can accomplish this with some kind of SQL select function in GRASS; however, I do not understand how the data is organized enough to do this.  When I import the SDTS file, I get 10 database links.  Most of the databases look 'empty', so I imagine there is just one file I need to be working with.

Do any of you have any experience in working with Hypsographic data?  Do you know how I can simply delete the entries that do not correspond to the lines I want?


texts: Highlighting area in Google Maps on mouse hover

texts:
I want to implement a feature in my app (gmap powered) to highlight the area on mouse over. Much like the functionality in wikimapia.  Check out the pic
or url

But i could't find any info on how to do this. I believe we need to have the shape info of every location to highlight the are, but i dont know how to do.

Can someone help me out?


texts: Increasing raster extent by filling extent with no data using ArcGIS Desktop without Spatial Analyst?

texts:
I am using feature to raster to rasterize water polygons. It only makes a raster the extent of the features in the polygons, whereas I would like it to be the extent of another existing .tif image.

Global Mapper can export a raster with a larger extent than the input and fills the extra space with no data, but I don't want to integrate 2 separate processes for a large number of files, I want to work in ArcGIS.

I am using ArcGIS Desktop 9.3 but I do not have the Spatial Analyst extension.

DOES NOT WORK IN ArcGIS 10.0, it won't respect the new envelope.


texts: When trying to make changes to a PostgreSQL view I get the following errors. Can someone advise how to get this working? I'm not familiar with triggers.


texts: Model works for single layer but errors for batch process

texts:
I have a pretty simple model in QGIS that I am trying to use as a batch process.

The model takes a line shapefile input, creates a temporary buffer, extracts by location using that created buffer, and then merges the extracted location shapefile with originally inputted line file.

When I execute it on a single layer it works fine and outputs as desired, however, when I try to use it as a batch process it errors out.

The Python error that kicks up:

  Traceback (most recent call last):

  File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python/plugins\processing\gui\BatchAlgorithmDialog.py", line 88, in runAlgorithm parameters[param.name()] = wrapper.parameterValue()AttributeError: 'NoneType' object has no attribute 'parameterValue'

And the model:


texts: how to snap two line vertex in two different layers used in QGIS

texts:
I want to snap between two poly-lines from 2 different layers (Layer1 in Brown and Layer 2(postGIS layer) in Pink color- ), but my snap option don't work. I have changed the tolerance to different levels but still not getting the results.

I don't get the magenta + sign to appear near vertices for snapping as many say it online.

Which tool used for snapping? I know only the "Node tool" which is not working .
I have set both layers editable (toggle editing set: ON)
Is there any better plugins there(I have tried many available ones.)?

Please give me suggestions since I am a beginner in QGIS.


texts: Why the layer can be loaded to Google Maps not OpenStreetMap using MapServer?

texts:
I am using OpenLayers WMS to load a layer from the PostGIS Raster Table, the JavaScript on OpenLayers are as following:

function init(){
    var map = new OpenLayers.Map('map');
    var imagery = new OpenLayers.Layer.WMS(
        "Global Imagery",
        "http://maps.opengeo.org/geowebcache/service/wms",
        {
            layers: "bluemarble",
        }
     );
    map.addLayer(imagery);
    localWMS = new OpenLayers.Layer.WMS(
        "Local WMS layers",
        "http://127.0.0.1/cgi-bin/mapserv.exe?MAP=C:/Users/gyu1/Desktop/Newfolder/mapserver.map",
        {
            layers: 'polygon',
            transparent: true
        }, {
            isBaseLayer: false
        }
    );

    map.addLayer(localWMS);
    map.setCenter(
        new OpenLayers.LonLat(1.147, 52.472).transform(
            new OpenLayers.Projection("EPSG:4326"),
            map.getProjectionObject()
        ), 5
    );
    map.addControl( new OpenLayers.Control.LayerSwitcher());
}

The map file linking the Data and styling are as following:

MAP
    NAME "polygon"
    WEB
        METADATA
            "wms_title" "MapServer"
            "wms_onlineresource" "http://127.0.0.1/cgi-bin/mapserv.exe?MAP=C:/Users/gyu1/Desktop/Newfolder/mapserver.map"
            "wms_srs" "EPSG:4326 EPSG:3857"
        END
    END
    PROJECTION
        "init=epsg:3857"
    END
    LAYER
        NAME "postgis"
        TYPE raster
        STATUS ON
        DATA "PG:host=localhost port=5432 dbname='fusiongeo' user='postgres' password='**' schema='public' table='london' mode='2'"
        PROCESSING "NODATA=-9999"
        PROCESSING "SCALE=AUTO"
        CLASS
            NAME "boring"
            EXPRESSION ([pixel] < 50)
            COLOR 250 100 0
        END
        CLASS
            NAME "mildly interesting"
            EXPRESSION ([pixel] > 50 AND [pixel] < 100)
            COLOR 100 0 0
        END
        CLASS
            NAME "very interesting"
            EXPRESSION ([pixel] >= 100)
            COLOR 0 50 0
        END
    END
END # end map

It works for the Google Base Map not the OpenStreetMap. I am quite new with the MapServer, OpenLayers and I don't know how to change the SRS to make it work.

The data in PostGIS are in EPSG:3857.


texts:
[I'm new to the site and my English is poor, so it is difficult for me to explain my problem. But there is no other person I can ask.]

My question is related to this exercise http://postgis.net/workshops/postgis-intro/joins.html, but it gives an error when I run it.

SELECT
subways.name AS subway_name,
neighborhoods.name AS neighborhood_name,
neighborhoods.boroname AS borough
FROM nyc_neighborhoods AS neighborhoods
JOIN nyc_subway_stations AS subways
ON ST_Contains(neighborhoods.geom, subways.geom)
WHERE subways.name = 'Broad St';

The error is Operation on mixed SRID geometries


texts: How to inverse the orientation of the holes coordinates of a polygon in geopandas

texts:
I have a polygon shapefile which I am reading with geopandas. The polygon contains holes. My purpose is to add a random point inside each hole.
In order to create the points inside the holes (interior polygon rings), the definition of their coordinates needs to be reversed (from counterclockwise to clockwise) otherwise the point falls outside of the polygon.
I get the coordinates of the interior rings using GeoSeries.interiors and what I get is

  LINEARRING (85002.811 446988.023, 85010.79399999999 446992.869, 85005.61900000001 447001.417, 84982.78599999999 446987.513, 84987.933 446978.99, 85001.89999999999 446987.47, 85002.811 446988.023)

How can I reverse the definition of the coordinates? Alternatively, I could get the x,y coordinates in a list and reverse it, but even in this case I cannot find a way to convert the geometry to list.


texts: Resources for learning to create mobile GIS Apps using Windows 8 Metro?

texts:
I have no experience with either mobile apps or Windows 8 Metro.

Are there any good resources out there that would be good for getting started from ground zero?


texts: Tilting osm map tile and other wms layer map tile like qgis map server tile via leaflet?

texts:
I have set up one osm map tile using mapnik and one wms map server using qgis map server. This both map server are able to render 2d map on web platform . I am now intrested to use these map server with custom looks. For now i am very much intrested to make the tile tilt with certain angle.

Can i make it using leaflet or by using other script?

I want to render my map tile with certain tilt angle like this.


texts: Thiessen polygons extend outside of extent of projection

texts:
I'm working on a geological project on the past position of continents (in ArcMap 10 Service Pack 5, build 4400). The tool I'm creating in Model Builder uses Create Thiessen Polygons and all input data are in Winkel-Tripel projection. The Create Thiessen Polygons tool however creates polygons that fall outside of the extent of the projection (see image), also when I specify the desired extent in the parameters of the CTP tool.

The result is that the Thiessen polygons (i.e. the shapefile) cannot be used by any other tools later on, because

  the added layer has an extent that is not consistent with the
  associated spatial reference information

I've already tried clipping the Thiessen polygons shapefile, but not even that works (the Clip tool does not accept the Thiessen polygons as input, because of the inconsistent spatial reference).

Any suggestions on how to deal with this?


texts: ee.Number.toFloat is not a function in Google Earth Engine

texts:
I am getting an error trying to parse a number to float in GEE, here is my code:

var getUHI = function(feature) {
  var urban = ee.Number(feature.get('LST_urb_day_CT'));
  var UHIday = ee.Number.toFloat(urban.subtract(ee.Number(feature.getNumber('LST_rur_day'))));
  return feature.set({UHIday: UHIday});
};

var UHIAdded = urb.map(getUHI);

currently in the featureCollection UHIAdded, 'LST_urb_day_CT' and 'LST_rur_day' are both type float, but after I parse them to numbers to calculate UHIday, UHIday becomes Number. I want to parse UHIday back to float but it tells me toFloat is not a function?


texts:
I make refactoring from Oracle to PostGIS and I found out that PostGIS doesn't support sdo_join.

How should I correct the statement below?

SELECT e.ID        AS EDGE1_ID,
       u.NODE_ID   AS NODE_ID,
       e.GEOM      AS EDGE1_GEOM,
       e.FTID      AS EDGE1_FTID,
       e.FID       AS EDGE1_FID,
       u.NODE_GEOM AS NODE_GEOM
FROM GRAPH_EDGE e,
     (SELECT j1.ID AS EDGE1_ID, j2.ID AS NODE_ID, j2.GEOM AS NODE_GEOM
      FROM TMP_EDGE j1,
           TMP_NODE j2,
          TABLE(sdo_join('TMP_EDGE', 'GEOM','TMP_NODE', 'GEOM', 'mask=EQUAL+CONTAINS')) x
   WHERE x.rowid1 = j1.ROWID
     AND x.rowid2 = j2.ROWID
     AND j1.ID != j2.ID) u
WHERE e.ID = u.EDGE1_ID


texts:
I have a few thousand TIF tiles and need to extract polygons of the actual data coverage from inside these TIFs (see below for image showing what the coverage looks like). The resulting geometry will be loaded into a Postgres table.

Due to the number of TIFs, and the fact that the data-set will be regularly updated, I need the process to be as efficient and automated as possible.

I have tried gdal_polygonize.py, but this gives me a result like this:

I would then have to 'dissolve' these polygons, or run ST_Union with PostGIS to get the result I want, which looks like this:

Dissolving is very slow, and feels unnecessary when there must surely be a way to skip the multiple polygon result of gdal_polygonize.

Does anyone have any ideas how I could do it in one step, or just make the process faster/more efficient?


texts: Opening NetCDF file in QGIS with corresponding coordinate system

texts:
I have a NetCDF file with 3 variables:

a 2D lat grid (degrees)
a 2D lon grid (degrees)
a 2D data grid

I am having troubles opening this file in QGIS with its corresponding latitude and longitude coordinates (see pictures below).
Format:
           netcdf4_classic
Dimensions:
           X12_251 = 240
           Y20_465 = 446
Variables:
    LON
           Size:       240x446
           Dimensions: X12_251,Y20_465
           Datatype:   single
           Attributes:
                       _FillValue = NaN
    LAT
           Size:       240x446
           Dimensions: X12_251,Y20_465
           Datatype:   single
           Attributes:
                       _FillValue = NaN
    SMB
           Size:       240x446
           Dimensions: X12_251,Y20_465
           Datatype:   single
           Attributes:
                       _FillValue = NaN

I am however able to open it in Matlab with the following script:
clear all
close all

ncdisp('nc_example.nc');
Z = ncread('nc_example.nc','SMB');
lat = ncread('nc_example.nc','LAT');
lon = ncread('nc_example.nc','LON');

geoshow(lat, lon, Z, 'DisplayType', 'texturemap')

I created a Google Drive with the .nc file here:
https://drive.google.com/file/d/1pbmmDjMaw0eWM1XIBk_px3piaAUtDi2f/view?usp=sharing


texts: Error installing resource from QGIS Resource Sharing Plugin

texts:
I have installed the QGIS Resource Sharing plugin in QGIS.

When I try to install the Mapbox's Maki Icons resource or OpenStreetMap style resource I get the following error message:

  Error: urlopen error Tunnel connection failed: 407 authenticationrequired

My network options in QGIS are currently set to Use proxy for web access with DefaultProxy set as the Proxy type.


texts: Using QGIS to prepare/execute fuzzy merge of two data sets with spatially coincident features?

texts:
I have two datasets that I need to prepare for merging: both are of a set of POI in the UK (places of worship) which have been generated by two different sources. There are around 50,000 points total that I am working with. I am aware that ArcGIS has "find identical" which can generate a combined table of geospatially coincident features to prep for merging, so I am wondering if there is a good way to do this in QGIS. I note that there is a solution here: Identifying duplicate features by comparing two layers in QGIS?, but this seems pretty dated already (QGIS 2.8.2 "Select by Location" dialog lacks "Precision" field). Additionally, I would like to add a secondary (and thus probably best to be very fuzzy) level of comparison on the name fields.

Sample csv snippet 1:

X,Y,UNIQUE_REFERENCE_NUMBER,NAME
64603,796126,17128583,Our Lady of the Waves and St John R C Church
65742,801811,17125921,St Brendan R C Church
66698,798282,68805407,Lady of Star of Sea
67070,803420,17125920,Barra Parish Church
70334,807612,17125924,St Vincent's R C Church
70696,803097,17125925,St Barr R C Church
72626,946839,20133561,Chapel
72687,870275,17125963,Kilmuir Parish Church
73286,829700,17176145,St Mary's Church
74062,868490,41252211,United Free Church
74496.6,821084.7,85602455,St Peter's Church
74911,868433,17125957,Free Presbyterian Church
75392,821388,41251293,Daliburgh Church of Scotland
75414,821411,68803787,Church of Scotland
75666,836413,17125941,Howmore Church
75816,816488,17125931,Our Lady of Sorrows Roman Catholic Church
75857.8,845895.8,41251641,St Michael's Church
76159,835684,17125940,St Joseph's R C Church
77416,841842,18117082,St Bride's R C Chapel
77881,852044,17125951,St Mary's R C Church
78641.3,812064.8,41249385,St Michael's Church

Sample csv snippet 2:

X,Y,NAME
67076.79711,803410.267,BARRA
72692.42456,870265.891,KILMUIR AND PAIBLE
75395.55481,821381.923,SOUTH UIST
75669.4918,836404.6101,SOUTH UIST

A test of this should match "Barra Parish Church" (group 1) with "Barra" (group 2) and should be fuzzy enough to match, "Howmore Church" to "South Uist" in the second instance in group 2.


texts: How to add data from a Google map as layer in ArcMap using C#?

texts:
I want to use C# to write a tool: add data from google map as a layer in ArcMap. I am  using ArcGIS 10.0.

Can you help me?


texts:
I have ArcInfo 10.2.1, Python 2.7 and FME 2013 SP4 installed on my Computer. I want to run a python script I exported from ModelBuilder (where it ran without issues) and imported in FME as shutdown script in another workbench model. While executing I get this message:

  The Esri Product license used by this machine does not include licenses for the 'ArcInfo' licensing level

I already tried the way specified here but this doesn't help with my problem. Any ideas what the error could be?


texts: How to style 2 columns with 900 combinations without resorting to rule based in QGIS?

texts:
So I have this point layer that has 2 columns that I need styled, one of which is owner_id and the other is size. Both are integer fields if it matters. The owner_id matches with an owner name in another table

The intended outcome is to have all 300 owners with names(instead of ids) to be listed by color, and for each owner a range of sizes to be styled by dot size. Essentially it would look like

**owner1name** green
size <5  green small
size 5-25 green medium
size >25 green big
**owner2name** blue
size <5 blue small
size 5-25 blue medium
size >25 blue big

Now I can do categorized style on the owner, switch to rule based and then add ranges but there are 2 problems with this approach

1 - Resulting ranges all take the same color instead of random ones

2 - If people toggle an owner on but untoggle one of the sizes, using the selection tool will also pick the untoggled combination

3 - Owners are listed by id, not name

What is the best way of tackling this?


texts: Value Map Attributes to Multiple Files in QGIS

texts:
I am hoping someone can share a way that I can some what automatically (none manual) associate a pre-defined CSV file for use with multiple shapefiles as an attribute Value Map.

I have a number (in the hundreds) of shapefiles that have the same attributes (Type1 and Type2).  Both Type1 and Type2 are numeric codes and I will like to apply text descriptions to them.  I have generated a CSV file that has the value and description.  What I am wondering is, is there a way to apply this CSV to all my shapefiles automatically, or will I have to open each file and associate the CSV file through the Properties dialog.

If there is an alternative method that would provide the same results I would be willing to give it a try.

Thanks...


texts: Extracting only roads through Overpass API query

texts:
I am using overpass API and querying the api for ways,as a result it gives points, whereas I need only lines, multiline strings and polygons. But while using the following code, I get a geojson file and when I load it into Qgis it gives me points and lines, when I switch off the points layer the lines are not loaded onto Qgis. I require only the roads(for which lines,multilines and multipolygons is required) and not Points. Kindly help as to how to modify the query to get only multilinestrings,lines and multipolygons and not Points.
import overpass
import geojson

api = overpass.API()

overpass_url="""

(
  way["highway"="motorway"](20.8,77.7,21.2,78);
  way["highway"="trunk"](20.8,77.7,21.2,78);
  way["highway"="primary"](20.8,77.7,21.2,78);
);
out body;
>;
out skel qt;

"""

res=api.get(overpass_url)

with open("test.geojson",mode="w") as f:
  geojson.dump(res,f)


texts: QGIS 3.0 Issue with rotation using Rectangles, ovals, diamonds (variable) processing tool

texts:
I'm trying to use the "Rectangles, ovals, diamonds (variable)" processing tool, using QGIS 3.0. One of the fields in my attribute table is a number going from 0 to 359 which is a reference to the "heading" direction of each point (meaning degrees..).

When I try to use this column as the Rotation field in the tool I get an error:

Traceback (most recent call last):
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\RectanglesOvalsDiamondsVariable.py", line 119, in processAlgorithm
self.rectangles(sink, source, width, height, rotation, feedback)
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\algs\qgis\RectanglesOvalsDiamondsVariable.py", line 145, in rectangles
feedback.pushInfo(QCoreApplication.translate('RectanglesOvalsDiamondsVariable', 'Feature {} has empty '
NameError: name 'QCoreApplication' is not defined

Execution failed after 0.01 seconds

It will accept various other fields as the Rotation field, most being random numbers. I've tried changing the type of my field from double to Integer and Int64 using the Refactoring Fields tool, but still the same issue.


texts: Viewer not displaying the data in my standalone PyQGIS Application

texts:
I tried creating a standalone external application of QGIS using python and I followed the tutorial mentioned in the link http://download.osgeo.org/qgis/doc/workshops/foss4g2007_qgis0.9_workshop_en.pdf

After running the python file the viewer opens and allows the user to browse for Shapefiles. The problem is that the Shapefile cannot be viewed in the viewer. It raises an error "Invalid Shapefile" though I added a correct Shapefile everytime I loaded any shapefile.

My code as well as the error is attached. Please tell me where I am wrong.

def addLayer(self):
    file = QFileDialog.getOpenFileName(self,
               "Open Shapefile", ".", "Shapefiles (*.shp)")
    fileInfo = QFileInfo(file)

    # Add the layer
    layer = QgsVectorLayer(file, fileInfo.fileName(), "ogr")

    if not layer.isValid():
        raise IOError("Invalid Shapefile")

    # Change the color of the layer to gray
    symbols = layer.renderer().symbols()
    symbol = symbols[0]
    symbol.setFillColor(QColor.fromRgb(192, 192, 192))

    # Add layer to the registry
    QgsMapLayerRegistry.instance().addMapLayer(layer)

    # Set extent to the extent of our layer
    self.canvas.setExtent(layer.extent())

    # Set up the map canvas layer set
    c1 = QgsMapCanvasLayer(layer)
    layers = [c1]
    self.canvas.setLayerSet(layers)

And the error window is


texts: Identify the geometry type of layer with QGIS expressions

texts:
The question
Using QGIS expression, it can occur that you want to identify the geometry-type of a layer (point, line or polygon? Multiple or single part?). This can be the case to write an expression that works for different kind of geometry types, using an if-clause - pseudocode: if point, do A, if line, do B, if polygon, do C.
See here for an example (go to the expression at B. Get list of all the angles of a line or polygon).
What I tried
I looked if there is a variable of the kind @layer_type, but there is not. There are variables for @layer_name or @layer_crs, but not for the geometry type of the layer (see screenshot).
As in the example above, a workaround is to find out if a geometry is a polygon or not to use if (num_rings( $geometry) is null, 1,0 ): returns 0 for polygons, 1 for lines and points.


texts: Printing name of each layer in mxd using ArcPy?

texts:
My ArcMap verison is 10.2.

My goal is to print the name of each layer in a mxd using Arcpy.

My script is below. There is no error, but it prints the first layer multiple times instead of all the layers in the "Main Map" data frame.

Does anyone have a solution?

import arcpy
mxd = arcpy.mapping.MapDocument(r"O:\Solar Databases\MO_projects.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Main Map")[0]
for lyr in df:
    x = 0
    name = arcpy.mapping.ListLayers(mxd, "", df)[x].name
    print name
    x = x+1
del mxd