texts: Is there a way to call the "Select By Location" dialog box and gain access to the parameters

texts: In an application that I am developing using VB.Net with ArcObjects 9.3, I am trying to display the "Select By Location" dialog box in my code.  I know that I can call it directly using the UID through the ICommandItem interface, but I also need access to the parameters set by the dialog box (e.g. the selection layer, the selection method, and so on-- something that I don't think is possible with the UID method) .  Is there a  way to do this somehow that I'm missing, similar to the way that you can call the "Select By Attribute" dialog box using the IQueryAttributes interface?


texts: Counting number of pixels from MODIS Burned Area Monthly Global 500m within shapefile boundary in Google Earth Engine

texts:
I have an ImageCollection derived from the MCD64A1.006 MODIS Burned Area Monthly Global 500m dataset available in Google Earth Engine. I've clipped it to the boundary of a national forest in California:

My goal is to count the number of red pixels within the boundaries of that national forest shapefile and then the total number of pixels within the shapefile to calculate the proportion of area burned. I've tried using code from a previous project, excerpted below, that would allow me to count each burned pixel with a system:index property, which they all have:

var burnStat = burnedSGSB.reduceRegions({
  collection: SGSB,
  reducer: ee.Reducer.countEvery().group({
    groupName: 'system:index',
  }),
  scale: 500
});

I suspect this has something to do with there being many images associated with that given pixel, but: What is the correct reducer needed to derive my target information, and what is preventing the above code from working?

Link to the full code is here:
https://code.earthengine.google.com/52917210ff19b3f6c300ce7d295e2261


texts: VBScript Syntax Error in ArcGIS labeling?

texts:
I have been using a script on older ArcGIS to convert depths in engineering format to the Hydrographic format seen on charts.  On trying to use it I get a syntax error in Line 14.  No matter how I look I cannot see/find the error.

Function FindLabel ([DEPTH])
  dim strSounding
  dim lngTrim
  dim lngDryHeight
  dim strSub

  lngDryHeight = instr([DEPTH],"-")
  lngTrim = instr([DEPTH],".")
  if lngTrim > 0 then
     strSub = right([DEPTH],1)
     strSounding = abs(left([DEPTH],lngTrim))
     if lngDryHeight > 0 then
        strSounding = "<UND>" & strSounding & "</UND>"
     end if
     strSounding = "<FNT name = ""Arial"" size=""7"">" & strSounding & "<SUB>" & strSub & "</SUB>" & "</FNT>"
  else
     if lngDryHeight > 0 then
        strSounding = "<UND>" & abs([DEPTH])  & "</UND>"
     else
        strSounding = abs([DEPTH])
     end if
     strSounding = "<FNT name = ""Arial"" size=""7"">" & strSounding & "</FNT>"
  end if

  FindLabel = "<ITA>" & strSounding & "</ITA>"
End Function


texts: Automatically assigning bins for vector graduated classes within a given range (QGIS 3)

texts:
Is there a way to assign graduated bins automatically for vectors within a given range in QGIS3?

I am trying to produce a series of maps which share a same legend, but do not necessarily exhibit all the same classes within each of their ranges:

For this first map, I assigned bins manually so that the entire range would be covered, but this is a time consuming process that I would have to repeat for each class and for each set of maps.

QGIS seems to automatically bin ranges from the smallest to the largest value present within a particular attribute, but I'd like to know if there is a way to override that so that I can provide a lower and upper value and automatically bin them.


texts:
I have two rasters, one for road network and another one for the DEM of the country.

I have to overlay the road network raster on top of DEM raster so as to get the elevation of roads. But my rasters have very different values and warp and resample didn't help.

This is my road network raster

And this is the DEM raster

Any ideas how I can accomplish my goal?

Both the rasters have same projection

'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]'

GDAL INFO of rasters

Here's the gdalinfo on the rasters :

Road Network Raster:

Driver: GTiff/GeoTIFF

Files: Samoa_reprojected.tif

Size is 1379, 607

Coordinate System is:

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]

Origin = (-172.802650000000000,-13.440900899999999)

Pixel Size = (0.001000000000000,-0.001000000000000)

Metadata:

  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=WhiteboxTools
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72

Image Structure Metadata:

  INTERLEAVE=BAND

Corner Coordinates:

Upper Left  (-172.8026500, -13.4409009) (172d48' 9.54"W, 13d26'27.24"S)
Lower Left  (-172.8026500, -14.0479009) (172d48' 9.54"W, 14d 2'52.44"S)
Upper Right (-171.4236500, -13.4409009) (171d25'25.14"W, 13d26'27.24"S)
Lower Right (-171.4236500, -14.0479009) (171d25'25.14"W, 14d 2'52.44"S)
Center      (-172.1131500, -13.7444009) (172d 6'47.34"W, 13d44'39.84"S)

Band 1 Block=1379x1 Type=Float32, ColorInterp=Gray

  NoData Value=-32768

DEM Raster:

Driver: GTiff/GeoTIFF

Files: SamoaElev_reprojected.tif

Size is 1204, 439

Coordinate System is:

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]

Origin = (-173.832990217792712,-13.115807779999948)

Pixel Size = (0.002301515790457,-0.002301515790457)

Metadata:

  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=WhiteboxTools
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72

Image Structure Metadata:

  INTERLEAVE=BAND

Corner Coordinates:

Upper Left  (-173.8329902, -13.1158078) (173d49'58.76"W, 13d 6'56.91"S)
Lower Left  (-173.8329902, -14.1261732) (173d49'58.76"W, 14d 7'34.22"S)
Upper Right (-171.0619652, -13.1158078) (171d 3'43.07"W, 13d 6'56.91"S)
Lower Right (-171.0619652, -14.1261732) (171d 3'43.07"W, 14d 7'34.22"S)
Center      (-172.4474777, -13.6209905) (172d26'50.92"W, 13d37'15.57"S)

Band 1 Block=1204x1 Type=Float32, ColorInterp=Gray

  NoData Value=-32768


texts: QChainage only works on first part of merged line

texts:
I'm trying to merge a number of lines into a single line for chainage purposes. Each individual line is created using OpenRouteService, then merged into a longer continuous shape using Edit > Merge Selected Features. However, when I then run QChainage on the resulting merged line, the chainage stops at the end of the first line, even though I'm using the merged shape.

See the screenshot below for an example. I've merged two lines, one north of Boyle, the other starting in Boyle and heading south. The chainage only appears on the northern section of the merged line.

Any ideas?


texts: Cannot set boolean operator within ArcMap ModelBuilder as True

texts:
I'm trying to set up my own tool to automate processing ASCII data (3D laser scanning data exported from the software CloudCompare) and have thus far got everything running pretty much as I would like it to. However, part of my workflow includes clipping an interpolated Raster using a shapefile [Clip (Data Management)].

When processing manually, I select "use input features for clipping geometry", however when I try to set this as an input variable in ModelBuilder, every time I click OK (after selecting the option or setting a boolean variable to True) it immediately reverts to 'False' (I can tell by hovering the mouse over the input).

Any ideas how to force the variable to stay as "True"?


texts:
I have 10 points that represent the centres of countries. I want to show with arrows the flow (some kind of trade) from one country to another.

The way that my table of attributes is constructed is the following:

From country      X Y       To country    X_to  Y_to            Value
1                 X Y       2             X_to  Y_to            10

I want the arrows to be proportional to the value and not to be straight lines but curved as it may be the case that 2 countries exchange trade, i.e. country 1 sends to country 2 and country 2 sends to country 1.

ARCGIS 10.2


texts: Debugging Visual Basic Code Block in calculate field using ArcGIS Desktop?

texts:
From time to time, I have to use Visual Basic in Calculate Field.

As the code get bigger, the more complicated is to debug it. And ArcGIS (9.3) does not help as the error messages are, in the best cases, something like:

"A field name was not found or there were unbalanced quotation marks.
ERROR 999999: Error executing function.
Failed to execute (Calculate Field)."

It's not very informative, has an ambiguous reference to the type of error and does not say the line of the code where the error occurred.

Is there a better way to debug the code in ArcGIS?

Does anyone have strategies to share, or a parallel software or site to do it?


texts: R -raster: Define raster extent from centre and rotation

texts:
I am quite new to R spatial transformation. My aim is to write a GeoTiff to be correctly geolocated in GIS software (e.g. ArcGIS).
I have a raster with known coordinates of the centre (proj: UTM WGS84  Zone 30 N; epsg: 32630). To my understanding (but please correct me if not), to display the raster in the correct position I need to define its extent (I know the pixel resolution). My problem is that the raster was originally captured from a UAVs platform so its correct position depend on the yaw angle of UAVs and gimbal at the time of shot.
How can I define correctly raster extension knowing coordinates of the raster centre and its yaw angle?


texts: Autoscrolling in ArcGIS Geoprocessing Results window?

texts: When checking the results window for messages while doing a lot of background processing it's rather annoying to scroll to the current message. So i was wondering if there is a trick to directly scroll to the current message line or is there such a feature in newer (>10.2) versions of ArcGIS?


texts: KML for Labels in Cesium to Display Countries, States, Roads, etc

texts: Is there a free way to display countries, states, roads, etc. in Cesium? I imagine these could be displayed using a comprehensive list of ground overlays in KML. I have looked for free sources of KML or KMZ to achieve these purposes, but have come up with no leads. Do you know of any place where I can get a list of country boundaries for instance in a KML file? Am I barking up the wrong tree? Are these labels are usually packaged with map tiles?


texts:
Someone mistakenly saved a GPS track as waypoints and then sent them along to me in a .csv file. They want a line back (shapefile). What's the simplest method to convert this to a line? Available tools are Arcmap, gdal/ogr and qgis more or less in order of preference. I'd rather not install an additional tool; an online conversion service would be okay.

Latitude,Longitude,Date,Time,ampm,,,
60.71586,-135.07476,25/07/2010,9:26:15,PM,,,,,
60.71637,-135.07563,25/07/2010,9:26:12,PM,,,,,
60.71648,-135.07612,25/07/2010,9:26:11,PM,,,,,
60.71664,-135.07707,25/07/2010,9:26:09,PM,,,,,
60.71672,-135.07756,25/07/2010,9:26:08,PM,,,,,


texts:
I have one column in my attribute table with attributes "1", "2" or "3".

1 = normal route
2 = special route because of the event (every year)
3 = special route because of the event and construction site (or other) (only in one particular year)

It's all about the routes of public transportation during a marathon race.

Is it somehow possible to create an note/remark for these attributes? I'm leaving the company soon and others maybe wouldn't understand my logic..;-)


texts:
I have made some python script to run under Processing toolbox.
A lot of these scripts use Gdal library and runs smoothly under windows and Linux.
But in mac os x 10.12.6 using qgis 2.18 the Gdal library is not exposed to the python console. I mean if i try

from os import system

system('ls')

I got 0, which means that it is ok
But if I type

system('gdal_translate')

I got 32512 which means command not found.
Of course if I run python from a terminal it is ok.


texts: Is there any reason not to upgrade to Windows 8.1?

texts: I have the option to upgrade to 8.1, but before I do I want to see if anyone has had any compatibility issues. I mainly run ESRI, but occasionally use QGIS, and other open source products. Has anyone ran into any compatibility issues with 8.1 and any GIS software?


texts: Building political world wall map for print?

texts:
I'm trying to create a political wall map. I've settled on the Natural Earth projection and thus have had problems with quite a few tools that don't support the projection (Tilemill, QGIS, Maperitive etc).

Trying to find out what the best approach is. Currently I think my best bet is to import my reprojected raster into Inkscape along with a raw SVG of my desired labels and do the label placement manually.

I am trying to generate proper SVG output from Mapnik for boundary lines and the raw labels for later adjustment, but so far I have been unsuccessful due to the size of the raster (21426x11141px). It simply crashes so far.

Since this is a personal project, investing in ArcGIS, Illustrator / MAPublisher seems a bit overkill.

Is there an easy way to convert my shp files into SVG for easy import into Inkscape - or is there a better approach than I am taking currently?


texts: Calculating Satellite Azimuth from Landsat WRS data?

texts:
I found the below answer in regards to how one identifies the satellite azimuth from Landsat Images, but am not sure how you translate the path/row into azimuth angle. I have searched for such, but am struggling to find an answer, therefore I repost in the hope that someone might be able to advise...

Alternatively, you could utilize the WRS (Worldwide Reference System) information (metadata fields WRS_PATH and WRS_ROW; for landsat 4 onwards, it is sort of a grid indicating where the satellite is on its orbit). Have a look at the USGS page with has WRS specification documents and shapefile downloads, it is basically straightforward.
However, if you are comparing scenes from the same area, the satellite's azimuth angle is not subject to any change.
(Note that Landsat Metadata has changed a little bit as of August 2012.)


texts: How to style a point based on Actual Attribute Values?

texts:
I have published a point layer in Geoserver, which is the data from a GPS Survey.
In addition to the geometry, my layer, also contains an Attribute for Accuracy, which indicates the Accuracy in meters of the reading.

I have been trying to style the layer, such that a point is show as a point surrounded by a circle whose radius is equivalent to the Accuracy in map units. (Similar to how most mobile apps show the GPS points)

All of the Attribute based symbolization that I have seen is Rule based (i.e. Classification).

How do you Symbolize a point, such that its radius is equal to its accuracy?


texts: Faster to polygonize a raster with internal compression or without? Python/GDAL

texts:
I'm trying to polygonize a bunch of rasters over the weekend (crunch time) and was curious what the thinking is regarding speed of polygonizing rasters in GDAL.

Based on your knowledge, if the raster had internal LZW compression, would the process be faster, slower, or the same? My rasters are very large (55500x55500 pixels) and since that will take a while to polygonize, I'm interested in making it as fast as possible.

I'm going to run a test on a subset but wanted to see if anyone who is familiar with compression (I'm not) and GDAL could answer this question.


texts: Unable to find a specific GRASS add-on in QGIS

texts:
I would like to use the r.xtent add-on made for GRASS GIS (http://svn.osgeo.org/grass/grass-addons/grass6/raster/r.xtent/description.html) in QGIS using the GRASS plugin, but I can't seem to find it.

What could be the issue?


texts: SIGSEGV when loading vector file with QgsVectorLayer class, yet ogrinfo looks valid?

texts:
I have a program set up as part of a QT application written in C++ to overlay data on a KDE Marble map, and am using the C++ API from QGIS to achieve this.  I create the vector file as a GeoJSON by launching a Python program using the QProcess class.  When I try to load the file as a vector layer, I get a segmentation fault error.  The line of code is:
const QString layer_name = "My Data";
const QString provider = "ogr";
QgsVectorLayer *vl = new QgsVectorLayer(file_path,layer_name,provider);

and the fatal returned to me according to GDB is from the internal function call:
QgsOgrProvider::getFeatures(QgsFeatureRequest const&) const()

and that's as far as the backtrace I have goes. Since I use Python to make the conversion into the JSON file and QGIS uses GDAL as the OGR provider, I check from the command line (Fedora Linux, V28) and I get no errors:
INFO: Open of `gdf.geojson'
      using driver `GeoJSON' successful.

Layer name: gdf
Geometry: Point
Feature Count: 1
Extent: (-83.475586, 21.579536) - (-83.475586, 21.579536)
Layer SRS WKT:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Geometry Column = _ogr_geometry_
U: Real (0.0)
V: Real (0.0)
T: Real (0.0)
P: Real (0.0)
spd: Real (0.0)
lon: Real (0.0)
lat: Real (0.0)
Time: String (0.0)
OGRFeature(gdf):96
  U (Real) = -9.57712459564209
  V (Real) = 7.50930976867676
  T (Real) = 21.077033996582
  P (Real) = 903.073608398438
  spd (Real) = 12.4758520126343
  lon (Real) = -83.4755859375
  lat (Real) = 21.5795364379883
  Time (String) = 2020-10-01_12:00:00
  POINT (-83.4755859375 21.5795364379883)

I'm not sure if maybe one of the features is invalid and is throwing it off or not as it's a large file (probably 100,000 point geometries; working on a GRASS application to make them into polygons and fix self-intersections and the like to make it more portable and easier to debug).  I've tried using rasters and meshes but I get other issues come up (the raster/mesh is from a NetCDF that I create from NCL).  I have the needed header files included in the C++ files, and my main.cpp initiates a QgsApplication:
QgsApplication b(argc, argv, true);
QString qgis_pfx = "/home/my_user/qgis";
b.setPrefixPath(qgis_pfx,true);
b.initQgis();

If anyone has any suggestions or needs to see more example code let me know.  I am using QT 5.11.3 with QGIS 3.8 and the Python version is 3.6.  Unfortunately most of the newer versions I have to build from source; upgrading to a newer Fedora OS is out of my hands.
EDIT: also I should add that when loading the file using the QGIS GUI, it works:


texts: How to import esri Personal Geodatabase (.mdb) to PostGIS on Linux

texts:
We have esri Personal Geodatabases that we need to import into PostGIS.

Ultimately, we'd like to allow people to upload an MDB on a web server, and have it import the layers into PostGIS.

On Windows, with FW Tools, the following works fine:

ogr2ogr -f "PostgreSQL" PG:"host=10.10.10.10 user=someuser dbname=poi password=somepassword port=5432" S:\GISData\Test.mdb -a_srs EPSG:26986

However, I'd like to automate this, (and do it programatically) so [I think] it would be a LOT easier on Linux.


texts: Merge all shapefiles in a directory - PyQGIS

texts:
I'm trying to merge all shapefiles in a directory using the following Python script from the console in QGIS 3.4.

The script was sourced from the following question:

Merging all vector files within directory with QGIS?

gis.stackexchange.com/questions/327398/merging-all-vector-files-within-directory-with-qgis

import glob, os, processing
path_to_shp = "P:/GIS/Basemapping/Topographic Area Schema 1/"
os.chdir(path_to_shp)
filelist = []
for fname in glob.glob("*.shp"):
    uri = path_to_shp + fname
    filelist.append(uri)

fileliststring = ';'.join(filelist)
processing.run("qgis:mergevectorlayers", fileliststring, "P:/GIS/Basemapping/Topographic Area Schema 1/merged.shp")

Running this results in:

TypeError: QgsProcessingAlgorithm.checkParameterValues(): argument 1 has unexpected type 'str'

The first time I ran into this I checked the input parameters for the algorithm using:

processing.algorithmHelp("qgis:mergevectorlayers")

Based on the algorithmHelp I assumed I should be passing a list of filepaths as input:

list[str]: list of layer sources

Running the following modification of the original script:

import glob, os, processing
path_to_shp = "P:/GIS/Basemapping/Topographic Area Schema 1/"
os.chdir(path_to_shp)
filelist = []
for fname in glob.glob("*.shp"):
    uri = path_to_shp + fname
    filelist.append(uri)

processing.run("qgis:mergevectorlayers", filelist, "P:/GIS/Basemapping/Topographic Area Schema 1/merged.shp")

results in the following error:

TypeError: QgsProcessingAlgorithm.checkParameterValues(): argument 1 has unexpected type 'list'

I feel I must be missing something obvious here. I've tried passing a list and a string of filepaths but neither seems to work. Can anyone help?


texts:
I have geographic data for the road network.

I changed my data from a geographic coordinate system (blue line) to a projected coordinate system (red line), it worked perfectly for all of my data, except for one year.

As you can see from the second image, a horizontal red line (road) appears from the centre south of the island and stretches off the land.

To try and identify the road, I used the select features tool. It only identified the red line on the land, when I select it off the land, nothing appears (like the red line is not there at all)

Do you know why this would happen? Or, how I can fix it? If I remove the red line, then I will be removing part of my data.

I would treat it as measurement error, but the first map appears does not have this error.


texts: Can't get dataset from raster map(.tiff) using GDAL on VS2012 C++

texts:
Can someone help me regarding this problem?

It seem the code can't read the data set of my raster map.It will throw this exception (popup to break or continue).

  Unhandled exception at 0x000000013F63298A in GdalTest.exe: 0xC0000005: Access violation reading location 0x0000000000000000.

// GdalTest.cpp : Defines the entry point for the console application.
//

#include "cpl_conv.h" // for CPLMalloc()
#include "stdafx.h"
#include "gdal.h"
#include <stdlib.h>
#include "iostream"
#include "gdal_priv.h"

using namespace std;

class GDALMajorObject;
class GDALDataset;
class GDALRasterBand;
class GDALDriver;
class GDALRasterAttributeTable;
class GDALProxyDataset;
class GDALProxyRasterBand;
class GDALAsyncReader;

int main(int argc, char** argv)
{
     GDALAllRegister();

    GDALDataset *poDataset;
    char* Raster = new char[256];
    double  adfGeoTransform[6];

    Raster = "@ C:\\Users\\User\\Desktop\\DSI REFERENCE\\Map Data\\Raster\\4257.tiff";

    poDataset = (GDALDataset *) GDALOpen(Raster, GA_ReadOnly )CPL_WARN_UNUSED_RESULT;
    if( poDataset == NULL )
    {
        cout<<"Damn it"<<endl;
    }

    cout<<( "Driver: %s/%s\n",
        poDataset->GetDriver()->GetDescription(),
        poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );

    cout<<( "Size is %dx%dx%d\n",
        poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
        poDataset->GetRasterCount() );

    if( poDataset->GetProjectionRef()  != NULL )
    cout<<( "Projection is `%s'\n", poDataset->GetProjectionRef() );

    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
        cout<<( "Origin = (%.6f,%.6f)\n",
            adfGeoTransform[0], adfGeoTransform[3] );
        cout<<( "Pixel Size = (%.6f,%.6f)\n",
            adfGeoTransform[1], adfGeoTransform[5] );
    }

    // make it 256 bytes for a string length of 255 plus null (\0) terminator

    // your code
    delete[] Raster;// release the memory VERY IMPORTANT

}

EDIT

It show this error when i try run the program.

Then I try to debug the code line by line. It got error reading character on string for this code line

char* Raster = new char[256];

when i about to debug another line it will jump to newaop.cpp tab.

// newaop -- operator new[](size_t) REPLACEABLE
#include <new>

void *__CRTDECL operator new[](size_t count) _THROW1(std::bad_alloc)
    {   // try to allocate count bytes for an array
    return (operator new(count));
    }

then i keep debug it will go debug to tab dbgdel.cpp until it stop at dbgheap.c then throw this

  memset.asm not found


texts: Proj.4/GDAL/QGIS Transformation between CRSs that are defined the same

texts:
I'm helping to ensure that open source software can handle Australia's new datum appropriately, see ICSM's website for details on the GDA2020 project.

Now, QGIS has already got the definitions of GDA2020 included, via GDAL, I understand.

An example GDA2020 coordinate reference system is this:

+proj=utm +zone=55 +south +ellps=GRS80 +units=m +no_defs

And if  you look at a GDA94 CRS, it's defined like this:

+proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

As you can see, these are very similar.

Now, the two CRSs are defined exactly the same, but, there's a shift in coordinates in GDA94 to GDA2020 of around 1.5 m to the north-east. (There's a grid shift file in NTv2 format that will soon be ready that will enable precise transformations, but that's not what this question is about.)

But, if you convert between GDA94 and GDA2020 now, using QGIS, there's no change in coordinates. It essentially just labels it differently.

Should there be a simple 7 parameter transformation implemented in Proj.4 or other open source tools that is the default transformation (albeit, imperfect) between GDA94 and GDA2020?

Or is it simply the case that the tools will always do no change?

How should this be handled?

(And I just want to note again that transforming using a grid is ideal, and that's handled in a few ways including this QGIS plugin.)


texts: What permissions does the CartoDB OAuth API provide?

texts: As a developer, I'm interested in finding out more about the CartoDB OAuth API. My sense is that the OAuth API permits the developer to make any SQL call on behalf of the user on any table, but I'm hoping it has a fine-grained (and as yet undocumented) permission system. My use case requires asking users permission to use some tables (but not all) and ideally it would let me get specify the type of access, i.e. read, write or read/write. Can I do that?


texts:
I have a point layer from an aerial photographic survey of a savanna, where each point contains information on the number of trees found within each photo. Some points can have 0 trees, while other points can have 1 or more trees.

Points representing photographs along an aerial survey. Each point contains the number of trees within each photograph

I generated a line layer from the points and split it into 1 km features. I wanted to sum the number of trees within each 1 km feature and then display the number of trees/km with a color (through quartiles, Jenks natural breaks, or some other classifier).

A line layer split into 1 km feature, with each feature containing a color indicating the number of trees based on the sum of trees from the points within each 1 km feature

I then wanted to find the average of five of these segments to get the average number of trees/km for 5 km features. Ideally this could be done to find the average number of trees/km for 5, 10, or even 100 km features.

Mean trees/km every 5 km

It seems like a very easy thing to do but I haven't figured out how to do it. I tried creating a raster file that had the sum of trees for 1 km cells, followed by converting that back to points and then converting that back to a raster with the average number of trees/km in 5 km cells. I'm not completely happy with that because the flight path isn't actually straight, so it would result in some raster cells represented by just one point and have extremely high or extremely low number of trees.

Is there a way to get the final result? I didn't want to use Kernal density since I wasn't trying to depict density over a unit area, but rather along a single linear transect. I've been using ArcMap 10.3.1 so far but I'd honestly be happy with using QGIS as well.


texts:
I am working on QGIS rasters.

What DEM file formats are supported by QGIS?

I know that QGIS can read the GDAL supported formats but not what the GDAL supported formats are.

I do not know whether QGIS supports all the DEMs supported by ArcGIS Desktop.


texts: Faster way to find nearest neighbors with LASlib?

texts:
I'm using LASlib to find the nearest neighbors of a subset of points in a point cloud.  I'm using brute force to find nearest neighbors, which is extremely slow (my code below).  Is there a better way, using the LASlib API to find nearest neighbors?

LASreadOpener lasreadopener;
lasreadopener.set_file_name("myLasFile.las");
LASreader* lasreader = lasreadopener.open();

float searchDistance = 1.0f;
I64 idx;
while (lasreader->read_point())
{
    //check to see if this point meets my conditions
    if (!conditions_met) // just pseudo code, i acutally check certain point attributes
    {
        continue;
    }

    idx = lasreader->p_count;
    double x = lasreader->point.get_x();
    double y = lasreader->point.get_y();

    // return to beginning
    lasreader->seek(0);

    // loop through all points
    while (lasreader->read_point())
    {
        if (lasreader->point.inside_circle(x, y, searchDistance))
        {
            // do something with this point
        }
    }

    //return to original loop point
    lasreader->seek(idx);
}
lasreader->close();
delete lasreader;


texts: Semi-transparent picture symbols in ArcObjects

texts:
When creating a new point layer with picture symbols in ArcMap 10.2.1 using ArcObjects, I have tried several ways to load semi-transparent symbol images from .png files, but to no avail, the transparency is always lost and appearing opaque.

Since IPictureMarkerSymbol.CreateMarkerSymbolFromFile doesn't handle .png, I think I will have to create the IPictureDisp myself. I know it should work somehow, since I can load the same semi-transparent .png file manually in the picture symbol in ArcMap, resulting in correct semi-transparency. And I know about the support for setting a transparency color, but it wouldn't handle semi-transparency.

My current solution, inspired from this post, goes something like this (except for the creation of the actual layer) but doesn't work either:

class IPictureDispConverter
{
  public static Guid iPictureDispGuid = typeof(stdole.IPictureDisp).GUID;

  /// Converts an Icon into a IPictureDisp
  public static IPictureDisp ConvertBitmapToIPictureDisp(Bitmap bitmap)
  {
    PICTDESC.Bitmap pictBit = new PICTDESC.Bitmap(bitmap);
    return OleCreatePictureIndirect(pictBit, ref iPictureDispGuid, true);
  }

  [DllImport("OleAut32.dll", EntryPoint = "OleCreatePictureIndirect", ExactSpelling = true, PreserveSig = false)]
  private static extern stdole.IPictureDisp OleCreatePictureIndirect([MarshalAs(UnmanagedType.AsAny)] object picdesc, ref Guid iid, bool fOwn);

  // WINFORMS COMMENT:
  // PICTDESC is a union in native, so we'll just
  // define different ones for the different types
  // the "unused" fields are there to make it the right
  // size, since the struct in native is as big as the biggest
  // union.
  private static class PICTDESC
  {
    //Picture Types
    public const short PICTYPE_UNINITIALIZED = -1;
    public const short PICTYPE_NONE = 0;
    public const short PICTYPE_BITMAP = 1;
    public const short PICTYPE_METAFILE = 2;
    public const short PICTYPE_ICON = 3;
    public const short PICTYPE_ENHMETAFILE = 4;

    [StructLayout(LayoutKind.Sequential)]
    public class Icon
    {
      internal int cbSizeOfStruct = Marshal.SizeOf(typeof(PICTDESC.Icon));
      internal int picType = PICTDESC.PICTYPE_ICON;
      internal IntPtr hicon = IntPtr.Zero;
      internal int unused1 = 0;
      internal int unused2 = 0;

      internal Icon(System.Drawing.Icon icon)
      {
        this.hicon = icon.ToBitmap().GetHicon();
      }
    }

    [StructLayout(LayoutKind.Sequential)]
    public class Bitmap
    {
      internal int cbSizeOfStruct = Marshal.SizeOf(typeof(PICTDESC.Bitmap));
      internal int picType = PICTDESC.PICTYPE_BITMAP;
      internal IntPtr hbitmap = IntPtr.Zero;
      internal IntPtr hpal = IntPtr.Zero;
      internal int unused = 0;

      internal Bitmap(System.Drawing.Bitmap bitmap)
      {
        this.hbitmap = bitmap.GetHbitmap();
      }
    }
  }
}

// FROM ANOTHER PART OF THE CODE:

IPictureMarkerSymbol markerSymbol = new PictureMarkerSymbolClass();
Bitmap bitmap = new Bitmap(pngFilename);
markerSymbol.Picture = IPictureDispConverter.ConvertBitmapToIPictureDisp(bitmap);

ISimpleRenderer renderer = new SimpleRendererClass();
((IGeoFeatureLayer)layer).Renderer = renderer as IFeatureRenderer;
renderer.Symbol = markerSymbol.Picture as ISymbol;

Here is the result with my code:

And here it is when I have pointed out the same .png file manually for the same layer in ArcMap afterwards:

Does anyone know how to do this properly?


texts: Renaming FCs in a Geodatabase with Python

texts:
I am new to Python but trying to rename hundreds of feature classes in a geodatabase. For instance, I have over 50 FCs named 'LAND_*'.

Can I with my Python script to run in ArcGIS Pro to rename the 'LAND' part of the file to 'INFRA'? Instead of 'LAND_A' you would rename to 'INFRA_A'?

I know how to do this individually but not how to do for all the feature classes in a geodatabase. I have looked through many examples but have not yet found one that works.

This is the code I have been using that did not work so far:

arcpy.env.workspace = "C:\\Temp\\Test.gdb"
for dataset in lstDatasets:
    lstFCs = arcpy.ListFeatureClasses("LAND_*", "", dataset)
    for fc in lstFCs:
        oldName = str(fc)
        newName = oldName.replace("LAND", "INFRA")
        arcpy.Rename_management(fc,newName)

I am getting the following error when I try to run the code:

"Traceback (most recent call last):
  File "<string>", line 1, in <module>
NameError: name 'lstDatasets' is not defined"

I was able to get to the exact solution I needed by using the following script:

arcpy.env.workspace = "C:\\Temp\\test_GDB.gdb"
lstDatasets = arcpy.ListDatasets()
for dataset in lstDatasets:
    print (dataset)
lstFCs = arcpy.ListFeatureClasses("LAND_*")
for fc in lstFCs:
    print (fc)
    oldName = str(fc)
    newName = oldName.replace("LAND_", "INFRA_")
    arcpy.Rename_management(fc,newName)


texts: Raster Images in Python to be ported to Swift in the future

texts:
I am working with a grey scale image from a DEM. No matter how much Googling or reading I do I can NOT find a solution for Getting elevation values from a raster images pixels. It should be as simple as referencing a legend and seeing that black equals X value and White equals Y value and interpolating the grey scale values in between. However this task seems impossible without libraries like GDAL or ArcGIS. I simple need to take an image in Python and grab a specific pixel then find the value within the pixel that from my understanding matches the grey scale value and a corresponding height.
I am looking to write this implementation in Python as I am stronger with the language and understand how to accomplish my goals using standard and additional libraries. I am writing a mobile application that mainly requires one thing... offline elevation data globally. I have my dataset that holds this information in raster images. I know the coordinates of every pixel. The only unknown is getting the elevation from each pixel. Each pixel is in grey scale, black being ZERO feet above sea level and white being the highest elevation. I am assuming that each intensity increase from one darker shade of grey to a lighter/more intense shade of grey is a 1 unit increase in elevation above sea level. I will correct for gravity later as I have completed this part. I simple need to know the implementation that takes the PNG file, then gives me the value of color for each pixel of which I can later change to coordinates and that's it. So I am avoiding libraries and will likely need to rewrite this code in Swift later when I decide to port this functionality over to the app. I am holding off on jumping directly into swift because I can complete a working example faster and which better understanding. I would not attempt to learn a new Math topic in a language I do not speak as I would be fighting the language and not grasping the topic.
How would I go about getting a single greyscale color value of a given pixel in an image in the Python language if I chose to avoid using common libraries and stuck to using only the Standard Python Library?


texts: export pymssql query to geodatabase table

texts:
I have thousands of databases in Microsoft SQL Server and I want the user to select the server and DB they want so I cannot use arcpy.ArcSDESQLExecute

I am trying to use pymssql to connect to the database and run a query. I understand how to save it as a CSV but I am struggling to understand how to save the results of the query to an in_memory table.

Here is my code to save it as a CSV.

from os import getenv
import pymssql
import pandas as pd
import arcpy
import csv
from decimal import Decimal

conServer = arcpy.GetParameterAsText(0) #r"server\server"
#conUser = getenv("put your username")
#conPassword = getenv("put your password")
condb = arcpy.GetParameterAsText(1) #"databaseName"

conn = pymssql.connect(server=conServer, database=condb)
cursor = conn.cursor()

query1 = "some sql query"

cursor.execute(query1)

with open(r"D:\GIS Data\output.csv", "w") as outfile:
    writer = csv.writer(outfile, quoting=csv.QUOTE_NONNUMERIC)
    writer.writerow([i[0] for i in cursor.description]) # write headers
    for row in cursor:
        writer.writerow(row)

cursor.close()
conn.close()

arcpy.SetParameter(3, outfile)


texts: Checking for overlaps in polygons using Oracle Spatial?

texts:
I want to get a list of all the polygons from table 2 which overlap with one polygon from table 1. Following is the query which I am trying to use which doesnt give any result.

 Select b.bg_id
 from table1 a,table2 b
 where a.store_id = '74073'
 AND SDO_Relate(a.geometry,b.geometry,'mask=anyinteract')='TRUE';

What would be wrong here?


texts: Is there a way to take an Extension License from another user in ArcMap/ArcCatalog?

texts: I'm using ArcMap 10.1 and need a specific extension that's activated on another user's machine.  Is there a script or other tool available to re-claim this license?


texts: Converting from Canada Atlas Lambert (ESRI:102002) to Lat/Long

texts:
I have a shape file with projected coordinates specified in Canada Atlas Lambert (NAD83) (ESRI:102002) and have been trying to convert the coordinates to lat/long (ideally in WGS84). However, whenever I do the conversion the map coordinates near the top of the map seem to disappear and the shapes continue straight to the edge of the map.

Notably, the map is of northern parts of Canada which seems to be the cause for the distortion.

Is there a way to properly convert the projection to coordinates? I've tried proj4 and QGIS with no luck.


texts:
I have a geoserver which I want to edit using c-sharp, I managed to retrieve information from it, but I haven't found a way to make any changes - I know that there is WFS-T but I haven't found a way to initialize it, I made the xml file (with my insert/delete/update), but I'm not sure how to send it to my geoserver, how can I send it using csharp?

URL = "http://localhost:8080/geoserver/topp"

XML file:

<?xml version="1.0" ?>
<wfs:Transaction
   version="2.0.0"
   service="WFS"
   xmlns:fes="http://www.opengis.net/fes/2.0"
   xmlns:wfs="http://www.opengis.net/wfs/2.0"
   xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
   xsi:schemaLocation="http://www.opengis.net/wfs/2.0
                       http://schemas.opengis.net/wfs/2.0.0/wfs.xsd">
   <wfs:Update typeName="topp:states">
      <wfs:Property>
         <wfs:ValueReference>STATE_NAME</wfs:ValueReference>
         <wfs:Value>Caliwhilenia</wfs:Value>
      </wfs:Property>
      <ogc:Filter>
         <ogc:PropertyName>STATE_NAME</ogc:PropertyName>
         <ogc:Literal>California</ogc:Literal>
      </ogc:Filter>
   </wfs:Update>
</wfs:Transaction>

request code

 WebRequest request = WebRequest.Create(url); //the geoserver url

    request.ContentType = "text/xml";
    request.Method = "POST";
    request.Credentials = new NetworkCredential("admin", "geoserver");

    byte[] buffer = File.ReadAllBytes(xmlrequest); //Makes a change to the state's name
    Stream reqstr = request.GetRequestStream();
    reqstr.Write(buffer, 0, buffer.Length);
    reqstr.Close();

    WebResponse response = request.GetResponse();

The geoserver is not changing the state's name

Update: now it connects to the server but I'm getting a "OperationParsingFailed" exception, altough I find that there is noting wrong with my xml file


texts: How to get all attributes from overlapping polygons with identical geometry into one single polygon in QGIS

texts:
I have one layer that has overlapping polygons with exact geometries. The most polygons that overlap are 5. For each polygon stack I would like to extract the values of a single field from EACH of the overlapping polygons. I don't want a summary value (mean, max, min etc), I want the actual values for every overlapping polygon. I could store these values in a shell polygon layer that had identical geometries if that helps. The dataset isn't massive (1484 features) but still too great for a manual approach.

I have searched and searched for various ways to do this. I can find ways to extract attributes from points to polygons and ways to summarise attribute values for overlapping polygons. But not to include the ACTUAL value of every overlapping polygon in an attribute table or something similar.

I have tried to union the layer with itself in QGIS. But this still results in overlapping polygons. I think it doesn't like that it has more than two overlapping polygons.

I have tried Join Attributes by Location, but this will only give me a summary value.

I have tried converting the polygons to rasters in case that path lead me down a solution. So I also have a folder with one raster for each polygon  in case that might help solve things. I created these using the v.to.rast.attr and selecting the input feature class to iterate through each polygon (little green arrow next to the input feature class field). But I couldn't get my head around a raster solution either.

Kazuhito provided a really easy solution. Here's how I ended up using it: I created the concatenated field using the expression given by Kazuhito, then deleted overlapping polygons and the other excess fields from the attribute table (no longer useful). I then saved the layer to .csv so that I could manipulate the concatenated column in excel. I extracted the comma separated into separate fields using the Data ribbon->text to columns tool. I then sorted the fields from smallest to highest value using a neat formula I found here https://www.pcmag.com/article2/0,2817,1215343,00.asp  I calculated a couple of extra columns that were useful to me and then imported the table back into QGIS as a .txt file and used a join (original layer properties ->joins) to add these attributes back into the layer. I'm sure there would be a way to do the excel manipulations in the QGIS field calculator. But for me this was easy enough.


texts: IWorkspaceEdit2Proxy equivalent in ArcGis Javascript api

texts:
Is there any equivalent for IWorkspaceEdit2Proxy in JS ArcGis API ? I've gone through the whole JS API documentation and haven't seen anything that would help to implement grouping edit operations.

I need to implement the same behaviour in Web application that is used in Desktop app. In ArcObjects SDK I can easily group edits and discard all of them later simple by using methods like IWorkspaceEdit2Proxy.stopEditing() IWorkspaceEdit.startEditing() and so on. However, it seems to me that ArcGIS Server manages edits applied on versioned dataset in his own way.

What can I do?


texts:
I try to create kml file from features.

My code

var feature=myVecLayer.getFeatureByFid('filedata.1');
                    var format;
                        function GetKMLFromFeatures(feature) {
                             format = new OpenLayers.Format.KML({
                                'maxDepth':10,
                                'extractStyles':true,
                                'internalProjection': new OpenLayers.Projection("EPSG:4326"),
                                'externalProjection': new OpenLayers.Projection("EPSG:4326")
                            });
                        alert(format);

                        alert(format);
                        return format.write(feature);
                        }

After i use this library to save file   https://github.com/dcneiner/Downloadify

                            Downloadify.create('downloadify',{
                            filename: function(){
                                                        alert("okey");
                                return  'kml_test.kml';
                            },
                            data: function(){
                                return GetKMLFromFeatures(features);
                            },
                            onComplete: function(){
                                alert('Your File Has Been Saved!');
                            },
                            onCancel: function(){
                                alert('You have cancelled the saving of this file.');
                            },
                            onError: function(){
                                alert('You must put something in the File Contents or there will be nothing to save!');
                            },
                            transparent: false,
                            swf: 'media/downloadify.swf',
                            downloadImage: 'images/download.png',
                            width: 100,
                            height: 30,
                            transparent: true,
                            append: false
                        });

But nothing heppens. If firebug i get error

TypeError: f.el is null

...tainer.id};if(f.options.enabled===false)f.enabled=false;if(f.options.transparent...

But more strange that i not get arert   alert("okey");

Whats wrong? Or if someone know better way to get kml file please give me advice.


texts: Would nearest point using Geodesic distance and nearest point using Haversine distance be the same point?

texts:
I have a point A and trying to find the nearest point to A in a list of points (B, C, D).

I could use knn with haversine metrics and get the nearest point like this:

knn = NearestNeighbors(n_neighbors=1, metric='haversine')

knn.fit(df['lat', 'lon'])

dist, idx = knn.kneighbors([(35.9157825, -79.0826045)])

However, I'm not sure if this point df.loc[idx] will always be the same point i'd get if I calculate distance using geodesic?

knn is very fast compared to having to calculate geodesic distance for all the points in my list. So I would love to use knn if the nearest point would always be the same.


texts: Can you reduce the number of contour lines?

texts: Using QGIS, is it possible to reduce the number of contour lines in a hypsography (topographic) layer? Mine has them every 5 meters, which is more than I need and it makes the map awfully busy.


texts: Interpolating Z at polyline vertices between known Zs

texts:
I have polyline strings with vertices of known elevation / z height and intermediate vertices with null/zero height. Known Zs are spread throughout the line.  I'd like to be able to populate the z of the intermediate vertices by linear interpolation between those known points.
I think this is functionality in Arc using IZ.InterpolateZsBetween?  Unable to verify due to no Arc access. Seeking functionality in QGIS or any alternative open package?


texts: Merging vector (building) elevations into DEM in QGIS?

texts:
I'm no expert in QGIS, I'm more into engineering.

I'm building a 2D hydraulic model and I would like to add the building footprints (elevated) to my DEM, so flow will go around structures in the model. I've actually done this before in QGIS but can't for the life of me get it to work again.

I've digitized buildings as polygons, added an arbitrary elevation to them, then tried a range of things to fuse that info with my DEM to no avail.


texts:
I am attempting to combine Lidar data that we capture using a scanner (which is real-world accurate) and USGS Lidar data. I am currently using Pdal and Entwine to download and reproject the section of USGS data I would like to use.
The problem I am foreseeing is that all the reprojections are still based on flat maps and don't take into account the curvature of the earth and how that would affect the height information. Is there a method in Pdal or other library/software that will correctly apply earth curvature to flat map derived Lidar data?
There might even be a way to apply a custom transformation to the data that would curve the position of the points based on the center of the bounding box?
Mostly I would ignore this distortion since we don't deal with large enough data sets that the loss of the correct horizon would cause a problem. But we want to potentially use mountains and other distant/large features accurately.


texts: zoomScale() doesn't work in QGIS 2.14 standalone c++ application

texts: Can't set scale for my canvas with mapCanvas.zoomScale(50000); Seems like canvas is zooming to 1:1 and for a few seconds I can't zoom in or zoom out with mouse wheel. Any suggestions?


texts: ArcGIS Server 10 - Feature Editing Services - Best Practices

texts:
What are some best practices for editing Multi-versioned sde data through a feature service?  What should I be looking into?

EDIT: seeing that when publishing a feature service the service is published but does not show up on the rest page..but I can reference it through it's full url??  If this observation would best serve as a separate question..please advise.


texts: Run QGIS layer Action on all selected features

texts:
QGIS is such an amazing tool with lots of useful features, so that it's easy to overlook some of them.

Just to be sure it's not already implemented in QGIS 2.18:

Is there any Attribute Table option to run a layer Action on all currently selected features?

I'm using a Processing script as workaround, but it would be nice to have this feature directly implemented in the Attribute Table.

##MyTools=group
##VectorLayer=vector
##action_name=string
from qgis.core import *
import qgis.utils
from qgis.gui import *

layer = processing.getObject(VectorLayer)
actionManager = layer.actions()
actions = actionManager.listActions()
if len(actions) > 0:
    cont = False
    i = -1
    for a in actions:
        i = i+1
        if a.name() == action_name:
            cont = True
            break

    if cont:
        selFeatures = layer.selectedFeatures()
        if len(selFeatures) > 0:
            for feature in selFeatures:
                actionManager.doActionFeature(0,feature)
        else:
            qgis.utils.iface.messageBar().pushMessage("Processing Error"," No features selected.",level= QgsMessageBar.WARNING,duration=5)
    else:
        qgis.utils.iface.messageBar().pushMessage("Processing Error"," Wrong layer action defined.",level= QgsMessageBar.WARNING,duration=5)
else:
    qgis.utils.iface.messageBar().pushMessage("Processing Error"," No layer action found.",level= QgsMessageBar.WARNING,duration=5)


texts: Extracting NDVI values from MODIS datasets in R

texts:
Firstly, I'll provide a bit of context behind my question. I've been reading a paper about bushfire modelling, and the paper uses the MODIS datasets (this is my first time learning about them). In particular, the paper uses fire detection data from the MCD14ML dataset, and NDVI values from the MYD13A3 (collection 6) dataset. I'm very new to MODIS datasets and HDF files. I've downloaded the 2003 MCD14ML dataset from https://earthdata.nasa.gov/earth-observation-data/near-real-time/firms/mcd14ml. Briefly, this dataset contains some attributes (including latitude, longitude, time) of detected "hotspots" in Australia, within the year of 2003. The dataset is in CSV form (so different attributes are in different columns). Below is an image of what I am seeing (in excel). So far, so good.

Unfortunately, the MYD13A3 (collection 6) dataset is not available in csv form (and the same applies for most other MODIS datasets). My goal is to add a column to the MCD14ML csv file which shows the relevant NDVI value for the region that was identified as a hotspot. So basically, something like this (I've made up the NDVI values).

I'm fairly confident that there's a way to download the MYD13A3 MODIS file in R, and that there's a way to extract NDVI values from it. Regarding downloading MODIS files in R, I tried to read the documentation for the R packages 'MODISTools' and 'MODIS'. Unfortunately, the documentations weren't very clear on how to use them to download MODIS files - so I searched on stack exchange, and found an answer that works at MODIS R Connection error . Based on the answer there, I've done the following.
hdf = getHdf("MYD13A3", collection = "006"
+              , extent="Australia"
+              , begin = "2003.01.01", end = "2003.01.03")

So what the above does is it downloads a bunch of HDF files to some specified folder on my computer. My understanding is that these HDF files contain the NDVI values of some geographic region. At this stage, what I want to do is to somehow extract the NDVI values from each HDF file, and put them into an R dataframe/CSV file. How can I do this? An additional question is whether there is a way (in R) to extract the NDVI values given longitude and latitude data.


texts: How do I uncheck a Toolbar button in ArcPy?

texts:
I have created a toolbar having 6 button to add/turn on different layers.

I want to have a functionality in my toolbar buttons, by which it will enable the user to see which layer is being added/turned on. Basically to turn the "checked" value as True.

Please see my code below.

class DG(object):
    """Implementation for DG_Imagery_addin.DG (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False

    def onClick(self):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        for lyr in arcpy.mapping.ListLayers(mxd):
            if lyr.longName in lyrlst:
                lyr.visible = False

        for lyr in arcpy.mapping.ListLayers(mxd):
            if lyr.longName == "MapProxy WMS\MapProxy WMS\DG":
                ispresnt = True
                islyr = lyr
                break
            else:
                ispresnt = False
        if ispresnt == True:
            islyr.visible = True
        else:
            df = arcpy.mapping.ListDataFrames(mxd)[0]
            imagery_lyr = arcpy.mapping.Layer(r"D:\Layer_Don'tDelete\DG.lyr")
            arcpy.mapping.AddLayer(df, imagery_lyr)
        arcpy.RefreshTOC()
        arcpy.RefreshActiveView()

class DG_CLEAR(object):
    """Implementation for DG_Imagery_addin.DG_CLEAR (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        for lyr in arcpy.mapping.ListLayers(mxd):
            if lyr.longName in lyrlst:
                lyr.visible = False

        for lyr in arcpy.mapping.ListLayers(mxd):
            if lyr.longName == "MapProxy WMS\MapProxy WMS\DG_CLEAR":
                ispresnt = True
                islyr = lyr
                break
            else:
                ispresnt = False
        if ispresnt == True:
            islyr.visible = True
        else:
            df = arcpy.mapping.ListDataFrames(mxd)[0]
            imagery_lyr = arcpy.mapping.Layer(r"D:\Layer_Don'tDelete\DG_CLEAR.lyr")
            arcpy.mapping.AddLayer(df, imagery_lyr)
        arcpy.RefreshTOC()
        arcpy.RefreshActiveView()

And 4 more buttons like this.

Now after clicking a button, how could I uncheck rest of the buttons and check my current button only.


texts: How to integrate Google Drive with web mapping capabilities?

texts:
I use a shared Google Drive spreadsheet to coordinate a large GIS mapping project.  The work involves processing thousands of DOQQ images.   Essentially, work is divided into study areas and then further into DOQQ boundaries that make up the study areas.

I am envisioning a dynamically updated web map that show which DOQQs have been processed within study areas.  In other words, as I update the shared spreadsheet, the web map should also update.  The attached image is a rough example of the type of data I wish to present.

Is this feat possible with any web mapping services (e.g. CartoDB)?  Which technologies should I pursue to accomplish this type of web mapping?


texts: How to get equidistant points from a linestring (geographical coordinates)


texts: What's happening with unusually large geopackage file size?

texts:
I have been working with a Geopackage in QGIS and R for a few weeks, storing ~ 12 layers of various geometry types.
Needing to email it to a colleague, I discovered the size of the file had grown to > 230MB
l *.gpkg
-rw-r--r--  1 user  staff   236M 24 Nov 12:37 AerialSurvey.gpkg

I was sure that there was nowhere near that amount of data in the system - so copied everything to a new geopackage:
ogr2ogr -f GPKG AerialSurvey-2.gpkg AerialSurvey.gpkg

The new copy is 1/30 the size!
l -h *.gpkg
-rw-r--r--  1 user  staff   236M 24 Nov 12:37 AerialSurvey.gpkg
-rw-r--r--  1 user  staff   6.7M 25 Nov 10:13 AerialSurvey-2.gpkg

My workflow involves analysis and re-saving of layers from time to time, from both QGIS and R.
Obviously this increase in file size is quick to fix with ogr2ogr, but why does this happen and how can it be prevented in future?


texts:
My questions is about heatmap plugin in QGIS 3.2.0

I have my data in CRTM05 projection, it´s a metric projection to Costa Rica based in WGS84 Datum.

When I try to run the Heatmap plugin, the running show this:

  "Processing algorithm… Algorithm 'Heatmap (Kernel Density Estimation)'
  starting… Input parameters: { 'DECAY' : 0, 'INPUT' :
  '/Users/robertovargasmasis/Desktop/Taller SIG 2018 CeU
  Neily/ejemplo_mapadecalor/puntos_crtm05.shp', 'KERNEL' : 0, 'OUTPUT' :
  '/var/folders/2z/fsqk6cy931d8r3c97qhc5gcr0000gn/T/processing_5bcc93e95ed54cd8b7e647b933ecf4e1/94fc27fbde38406eabee152a81d702aa/OUTPUT.tif',
  'OUTPUT_VALUE' : 0, 'PIXEL_SIZE' : 0.1, 'RADIUS' : 100, 'RADIUS_FIELD'
  : None, 'WEIGHT_FIELD' : None }

  Traceback (most recent call last): File
  "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/algs/qgis/Heatmap.py",
  line 215, in processAlgorithm self.tr('Could not create destination
  layer'))
  _core.QgsProcessingException: Could not create destination layer

  Execution failed after 0.03 seconds

  Loading resulting layers Algorithm 'Heatmap (Kernel Density
  Estimation)' finished"

Someone have had this problem?
How to fix it?


texts:
I'm trying to extract all data OSM within a GeoJSON polygon using osmium extract. I want the output to include only nodes and ways completely contained within the polygon. I tried using the --strategy simple, however this still includes ways that have at least one node within the polygon.

I checked quickly using the Overpass API and it has the same behaviour.

Is there a way to do this using osmium or osmosis? Or do I need to write my own program?


texts: R can't work with ecw while gdal/qgis does. Ubuntu 20.04 config

texts:
I'm an Ubuntu 20.04 user with all configured to use .ecw raster format. I can read these files in QGIS (3.10) and work with them without a problem. This is what I get in console:

$ gdalinfo --version
GDAL 3.0.4, released 2020/01/28

$ gdalinfo --formats | grep -i ecw
ECW -raster- (rw+): ERDAS Compressed Wavelets (SDK 5.4)
JP2ECW -raster,vector- (rw+v): ERDAS JPEG2000 (SDK 5.4)

The issue here is that I'm trying to perform some processes with this files within R but it seems that .ecw are not working and I don't know what to do.

> gdalbuildvrt(gdalfile = "PNOA_MA_OF_ETRS89_HU29_h50_0151.ecw", output.vrt = "dem.vrt")

ERROR 4: `PNOA_MA_OF_ETRS89_HU29_h50_0151.ecw' not recognized as a supported file format.
Warning 1: Can't open PNOA_MA_OF_ETRS89_HU29_h50_0151.ecw. Skipping it
NULL
Warning message:
In system(cmd, intern = TRUE) :
  comando ejecutado '"/usr/local/bin/gdalbuildvrt" "dem.vrt" "PNOA_MA_OF_ETRS89_HU29_h50_0151.ecw"' tiene estatus 1

What could I check? I prefer any solution that downgrading the gdal version to 2.4 (as is Ubuntu is quite complex to make it all work and that will brake other things for sure...).


texts: How to fix unusual initial DEM input values and derivative raster output values in SAGA GIS without rescaling?

texts:
I am trying to run elevation derivative tools on a provincial DEM provided by the Ontario Government website (Ministry of Natural Resources and Forestry) and I am getting strange values when bringing the data into SAGA GIS.

The values show up like normal in ArcMap but I get an unusually extremely low minimum value (that should be 0) as displayed in the image below:

As displayed, I am currently trying to create a Topographic Position Index (TPI) and the output values are also unusually large. Only when I adjust the min/max values can I get a proper visualization in the display windows, but the metadata in the "Description" tab does not change. On top of that, when I save the SAGA file and bring it into ArcMap I get more strange values. Only when I rescale using the "Rescale by function" (linear transformation function) tool in ArcMap can I get a normal visualization in ArcMap, similar to when I changed the raster's min/max values in SAGA and even then I am unsure if I input the parameters correctly into the ArcMap tool.

Would anyone also know what a typical standard TPI Index value range should be? And would anyone know the proper parameter inputs for a linear transformation function via the "Rescale by function" tool in ArcMap?


texts: Calling and showing photo from layer shown in composer?

texts:
What I'm trying to do is to see images from photo-layer that are shown in the composer.

There will be 6 images in each map. There will be a datatable from the point layer that contains the path and the name of the photo file.

What I want is that the first photo-square calls the first row from the table and shows the photo from the path. The same process for each of the squares.


texts: Force extent to ST_AsRaster() function in PostGIS

texts:
I am trying to create a raster from a PostGIS geometry, the issue is that wherever features are missing PostGIS tries to reduce the extent of the raster. But I do not want this behaviour.

Can I force extent in PostGIS raster, i am using the following code.

WITH q AS (
WITH a AS (
    SELECT
        ST_GeomFromText('POLYGON((1840136.2116530193 4984382.618753927, 1840136.2116530193 4989382.618753927, 1845136.2116530193 4989382.618753927, 1845136.2116530193 4984382.618753927, 1840136.2116530193 4984382.618753927))', 3857)
        AS geom
)
SELECT
    ST_Transform(a.geom, 4326)
    AS geom
FROM a
),

features AS (
    SELECT
    ST_Transform(ST_Intersection(t.wkb_geometry, q.geom), 3857)
    AS geom
    FROM
        public.gis_osm_roads_free_1 t,
        q
    WHERE
        ST_Intersects(t.wkb_geometry, q.geom)
        AND
        1=1
),

q1 AS (
    SELECT ST_AsRaster(
        ST_Collect(features.geom),
        100,
        100,
        '8BUI'::text,
        1,
        254,
        1840136.2116530193,
        4989382.618753927,
        0,
        0,
        false
    ) as ras
    FROM features
)

SELECT
ST_AsTIFF(
               q1.ras,
                'LZW'
            )
FROM q1


texts: Issues with spaces in file path in GRASS r.in.lidar function

texts:
Am using the r.in.lidar function in GRASS to manipulate multiple lidar files within GRASS.

For testing purposes i did everything with one file only by using the following line and this worked.

r.in.lidar -e -n -o --o output=Rst_Min_cf2  input="C:\Users\XXX\OneDrive - XXX\Projects\XXX\XXX_1K_Las.laz" method=min resolution=5 class_filter=2

However, i have more than 20 files that i would like to analyse together as file by file would be a bit cumbersome.

r.in.lidar allows the user to send a list of files to analyse together and i have successfully used that in the past. I tried to do the same in this case with the following line

r.in.lidar -e -n -o --overwrite output=Lidar file="C:\Users\XXX\OneDrive - XXX\Projects\XXX\AllLidarFiles.txt" method=min resolution=5 zrange=-1,100 class_filter=2

However, it did not work and i got

ERROR: Input file XXXX does not exist

After trying many thinks, i moved all my files to another folder making sure there were no spaces in the path

r.in.lidar -e -n -o -p --overwrite output=Lidar file=C:\Users\XXX\Desktop\XXX\AllLidarFiles.txt method=min resolution=5 zrange=-1,100 class_filter=2

and this did work, which made me think that the issue arises because of the spaces in the file path.

Does anyone know exactly why this problem occurs and how i could solve it without moving all my files to a folder whose path does not have any spaces at all?


texts: Should Identify Tool be used in a Web application?

texts: I was told having a Identify Tool is a bad practice in most web applications.  Can anyone explain why or point me to a resource that explains why having a Identify Tool is considered bad design in a web application?  If every feature you click on brings up a pop-up that could probably get annoying as well.  Seems like an identify tool could help prevent these potentially annoying pop-ups.


texts:
In QGIS, I use the python console to generate image files through the composer. In particular, I am creating a map with a legend. The legend is too big. It occupies about 1/8 of the image. I would like to make the legend smaller.

The code that adds the legend to the composition is:

c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)
cleg = QgsComposerLegend(c)
cleg.model().setLayerSet(mapRenderer.layerSet())
cleg.setItemPosition(0,45)
c.addItem(cleg)

A legend is created, but it is too big. In order to make it smaller, I have been reading the documentation of the QgsComposerLegend and QgsComposerItem classes, with no positive results. In particular, I have been trying to change the legend size by typing

cleg.setItemPosition(0,45,3,10)

but the size does not change.

I think I might have to make the font smaller first. If that is the case, I should use setStyleFont, which depends on a QgsComposerLegendStyle object which I don't know how to retrieve (I have tried to use the rstyle function, but I don't know how to use it).

In summary, could you please provide an example on how to:

Make the legend smaller
Make the legend font smaller?


texts:
I've tried to plot a map of magnitude earthquakes in Colombia with ggplot. After many efforts, I got it. But, I need to focus in specially area where the data of magnitudes is (4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5). When I use the same script for this data, I get this:

  "Error: Discrete value supplied to continuous scale".

I've used many solutions in blogs but it does not work.

This is my first map and the script

santos <- readOGR("Departamentos_WGS84.shp","Departamentos_WGS84")
map <- fortify(santos)
datos=read.table("coord.txt",header=T,dec=",",sep=";")
myData = data.frame(datos)

gg <- ggplot()
gg <- gg + geom_map(data=map, map=map,
                    aes(x=long, y=lat, map_id=id, group=group),
                    fill="#fffcfc", color="#0e0e0e", size=0.15) +
                    geom_point(data=myData,
                    aes(x=long, y=lat, size=Magnitud), color = "red",shape = 21) +
                    scale_size(range = c(0, 8)) +
                    coord_map() +
                    theme_bw() +
                    labs(x=NULL, y=NULL) +
                    theme(panel.grid=element_blank()) +
                    theme(panel.border=element_rect())
                    gg


texts: Using list of attributes to create layers?

texts:
I have a very long list of values (over 2000 rows) in the following format:

1: TAY100, AB60, EM10
2: CA40, 9, UT26, 24
3: YO56, YE

I have a layer with features each with a distinct attribute called VistaID that looks like this:

How can I create a layer with the name of the number on the left that includes the features with VistaIDs that correspond to the values in the list?

In other words, how can I create a layer called 1 that includes the features that have the VistaIDs of TAY100, AB60, and EM10, and then repeat this for every row?


texts: How to save space in the legend and write the layer name besides the symbol?

texts:
I'm starting using QGIS to make a project for university and I am with a problem on configuring legends on Print Composer.

The thing is I want the text of my legend to appear after the symbols/coloured rectangles and not above. My legend doesn't fit in the space because I have 2 lines for each layer when I could have just one. Is there any way I can solve this?


texts: How to Turn Standard Web Server into GIS Capable?

texts:
I have been pointed to help our IT dept to embeded GIS ability in our custom ERP system. The server that host the system has the following standard settings : Apache, PHP, MySQL and JavaScript. The problem is our IT guys not allowed me in "messing" with their system, so to suggesting installing mapserver or geoserver or other gis capable server is highly restricted.

Their current solution is converting our native mapinfo tab files into svg files. IMHO this solution will not answered the management needs.

I am thinking to use pure openlayers, another problem is the only way to read our tab files is by convert them into geojson.

The concept on their mind is all of GIS data is inserted into MySQL, I've heard about MySQL-GIS but I don't know how the plugin is supported today.

Basicly I don't do the job, I just give them the right direction, but with this limitation I can't figure out how to achieve this.

So to cleared out my question is, I have a standard Apache, PHP, MySQL and JavaScript system, I'm hoping my GIS data is read and write from or to MySQL. How can I turn that standard web system into GIS capable system without changing any setting on my server? OpenLayers are welcome.

Any suggestion are most welcome.

Regards,
Zery


texts:
I have two rasters (Hyrdological Soil Group (HSG) and Land Cover) and i'm trying to find the average HSG for each land cover type. The HSG raster contains the values 2,3 and 4 and there are 23 different land cover types. I've tried using zonal statistics using the mean statistics type, but this results in mean values ranging from 2.58 to 3.47.

I would like just single numbers (i.e 2,3 or 4) for each land cover type, is there a way around this?


texts: Could not locate LayerGroup in catalog GeoServer + OpenLayers?

texts:
I have a Layer Group in GeoServer which contains 4 regular layers (layers and group belonging to "TRAC" workspace). I'm trying to do the following on map clicking:

featureRequestGS = new WFS().writeGetFeature({
  srsName: 'EPSG:25830',
  featurePrefix: 'TRAC',
  featureTypes: ['groupLayer'],
  outputFormat: 'application/json',
  filter: containsFilter('geom', clickPoint, 'EPSG:25830'),
});

Once I make the request to GeoServer, I get that GeoServer couldn't locate groupLayer in catalog:

<?xml version="1.0" encoding="UTF-8"?><ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http://localhost:8080/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="InvalidParameterValue" locator="typeName">
<ows:ExceptionText>Could not locate {trac}groupLayer in catalog.</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

The funny thing is that on map initialization I set an WMS layer with the GeoServer layer group and it works fine:

this.groupLayer =
new ImageLayer({
  source: new ImageWMS({
    url: environment.baseUrl,
    params: {LAYERS: 'groupLayer'},
    projection: proyeccion
  })
});

this.map = new Map({
  target: 'map',
  layers: [this.groupLayer],
  view: new View({
    projection: proyeccion
  })
});

Any idea of what am I doing wrong here? I'm pretty stuck with this!

P.D: This is the GetFeature request body:

<?xml version="1.0" encoding="UTF-8"?>
<GetFeature xmlns="http://www.opengis.net/wfs" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" service="WFS" version="1.1.0" outputFormat="application/json" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
   <Query typeName="TRAC:groupLayer" srsName="EPSG:25830">
      <Filter xmlns="http://www.opengis.net/ogc">
         <Contains>
            <PropertyName>geom</PropertyName>
            <Point xmlns="http://www.opengis.net/gml" srsName="EPSG:25830">
               <pos srsDimension="2">60000 47000</pos>
            </Point>
         </Contains>
      </Filter>
   </Query>
</GetFeature>


texts:
I have been trying to figure out how to set the transparency of my labels and label buffers in PyQGIS for an application I am building. I am using the "setCustomProperty" class from the QGIS API. I can't find anything in the API documentation that allows for setting transparency. This is possible in the QGIS GUI via a slider, but seemingly impossible in PyQGIS. Here is a snippet of my labeling code:

    self.mylayer.setCustomProperty("labeling", "pal")
    self.mylayer.setCustomProperty("labeling/isExpression", True)
    self.mylayer.setCustomProperty("labeling/enabled", True)
    self.mylayer.setCustomProperty("labeling/bufferDraw", True)
    self.mylayer.setCustomProperty("labeling/bufferSize", "1")
    self.mylayer.setCustomProperty("labeling/bufferColor", "white")
    self.mylayer.setCustomProperty("labeling/fontFamily", "Arial")

Shouldn't there be a custom property setting for transparency? Like a bufferAlpha or a textAlpha setting or something?


texts: Rmarkdown build (index.html) in SublimeText 3 not displaying leaflet map, but `rmarkdown::run()` works

texts:
First Attempt with SublimeText 3 R Markdown Build
If I build with R Markdown in SublimeText 3, an index.html file is generated. If I then open this .html file in either Firefox, Chrome, Brave, etc, the map does not render. For example:
Result

Second Attempt with Terminal
However, If I switch to my Windows terminal and load R, then use rmarkdown::run(), the map generates in either browser of choice. For instance:
Result

Reprex Code: flexdashboard with shiny integration
---
title: "COVID-19 Tracker"
date: "`r Sys.Date()`"
fontsize: 9pt
fontfamily: Segoe UI
output: flexdashboard::flex_dashboard
runtime: shiny
---

suppressMessages( require(dplyr) )
suppressMessages( require(tidyr) )
suppressMessages( require(flexdashboard) )
suppressMessages( require(leaflet) )
suppressMessages( require(maps) )
suppressMessages( require(mapdata) )
suppressMessages( require(shiny) )
suppressMessages( require(htmltools) )
suppressMessages( require(htmlwidgets) )

cntry <- map("world2Hires",fill=F,plot=F)$names[1:100]
cntry <- cntry[order(cntry)]
getBounds <- function(vec=NULL)
{
    lng1 <- rep(0,length(vec))
    lat1 <- rep(0,length(vec))
    lng2 <- rep(0,length(vec))
    lat2 <- rep(0,length(vec))
    for (i in 1:length(vec))
    {
        lat1[i] <- map("world2Hires",fill=F,region=vec[i],plot=F)$range[3] %>% round(5)
        lng2[i] <- map("world2Hires",fill=F,region=vec[i],plot=F)$range[2] %>% round(5)
        lat2[i] <- map("world2Hires",fill=F,region=vec[i],plot=F)$range[4] %>% round(5)
        lng2[i] <- map("world2Hires",fill=F,region=vec[i],plot=F)$range[1] %>% round(5)
    }
    cbind(lng1,lat1,lng2,lat2)
}

location <- getBounds(vec=cntry)
df1 <- data.frame(cntry,location)

Column {.sidebar}
--------------------------------------------------

### Epicenter Region

# https://community.rstudio.com/t/rstudio-connect-interactive-markdown-with-shiny-error-path-for-html-dependency-not-provided/7263/11
selectInput("loc",
        label="Select Region:",
        selectize=FALSE,
        choices=df1[,1],
        selected=df1[19,1])

Column
--------------------------------------------------

### COVID-19 Epicenters

renderLeaflet({
    leaflet() %>% addTiles() %>%
    # add tile theme layer
    addProviderTiles(providers$CartoDB.Positron) %>%
    # add bounds
    fitBounds(lng1=map("world2Hires",fill=F,region=input$loc,plot=F)$range[1],
          lat1=map("world2Hires",fill=F,region=input$loc,plot=F)$range[3],
          lng2=map("world2Hires",fill=F,region=input$loc,plot=F)$range[2],
          lat2=map("world2Hires",fill=F,region=input$loc,plot=F)$range[4])
})


texts:
I need QGIS to display combo boxes when editing for certain fields in a table, very similar to the functionality provided by the Relations properties in the project settings.  It seems this is designed to simply have a value in one table refer to a list of possible values in another table, linked by foreign key.

However, my lookup table contains data for many different feature classes and more than one field between the two tables need to match.  Any given field of a feature may or may not have subtype/domain values (not sure the correct terminology), so it's not absolute that the columns in two different features from the same feature class (table) will contain combo boxes for the same columns.  It must be determined by simply querying the lookup table to see if that field has entries.

For any given feature, I need to query the lookup table for the following:

if the 'featureclass' column of the lookup table has a value for the current featureclass (table/layer name with the vector data)
if the 'code' column of the lookup table matches the 'code' column of the feature
if the value of the given feature column exists in the lookup table's "field_name" column

I've developed a plugin that does just that.  It loops through all of the fields of a selected feature and performs queries against my lookup table for each column.  A query looks something like this:

SELECT *
FROM lookup_table
WHERE code='ABC'
AND LOWER(field_name)=LOWER('def')
AND LOWER(featureclass)=LOWER('my_feature_class')

If the result of the query is empty, then I build a simple QLineEdit widget for editing the text directly.  If the result has values, I store them in a dictionary like this:

{ 1000: True, 1001: False, -999999: No 'Information' }

So where my feature has the value "1001" in the "def" column, I now build a QComboBox with the options "True", "False", and "No Information", and the combo box stores the integer codes as the data.  When a user makes a change and clicks "Save", I manually execute a SQL command to update that record with the new integer value.

Is this functionality possible natively in QGIS using relations?


texts:
I’m running ArcGIS Server 9.3.1 on a Windows Server 2008 virtual machine.  I noticed the other day that my C:\ drive was almost full.  I tracked the issue down to  this location “C:\Users\ArcGISSOC\AppData\Local\Temp” on the server.  There are some .mem files that are written everytime a service is accessed.  For the most part they are small, but every now and then one gets written that’s 2GB.  That is causing issues.  I’ve looked on the ESRI site and can’t find a resource to explain these files and what is being written to them.

Can anyone point me to a resource on these files that explains what they are, or have an explanation as to why I’m getting these large files?


texts: QGIS error message "raster layer data provider is invalid"

texts:
When I try to add a new WMS Layer from an URL into QGIS, it seems I can add the layer but when I try to open it I get the following error message: "invalid layer: raster layer data provider is invalid".

Any idea where the problem is? Is it the URL (which should be correct) or something with my QGIS. I'm using the QGIS 3.4.6 Version.

the WMS URL used was: "https://www.wms.nrw.de/gd/bk050?"

I have not used an other Provider because it is the only one I know for these data


texts: How to sort data on the closest distance basis in QGIS?

texts:
I'm preparing a map in QGIS 3.16.5 with the set of points and polygons in two layers. Unfortunately during a field work data were collected a bit randomly and objects of subsequent names are scattered on the map (by object I mean both: points & polygons). For instance, I have an object named 100 very close to object 1, while object 2 is located in a large distance. I want to rename objects to get some logical order of the names. I was going to export data to .csv, then change names in spreadsheet and import it again to a new layer.
Before that I need to sort data somehow using location and distance between objects as a rule. I would like to have object no. 1 in the first row, then, in the second row, the object that is closest to 1, in the third row object that is closest to this from the second row, in the fourth row object that is closest to this from the third row, and so on. In the end I would like to have numbers assigned respectively to each object.
Can it be done?


texts: How to update element text in arcpy.mapping?

texts:
I have a problem updating a map element text (named element1) in a map document using the following code:

a1 = '"' + ("element" + str(1)) + '"'
for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
if elm.name  == a1:
    elm.text = str("new text")
    elm.elementPositionX = 0.6046
    elm.elementPositionY = 6.4636

The code above does not update the element1 with the new text value.

It looks like the following part does not work:

if elm.name  == a1:

Because when using the below hard coded approach:

if elm.name  == "element1":

It worked, but I need to have it as variable.

What am I doing wrong here?


texts: Pixel coordinate to world coordinates using Python

texts:
I have developed a semantic segmentation method to map certain objects in aerial imagery. Throughout this project I download pictures (4800x4800 pxs) where I know the min_x, min_y, max_x and max_y in world coordinates, e.g. (6.212454957892032, 51.58908266914109, 6.219403242107967, 51.59339941284951). Resolution is 10 cm per pixel.

To then segment the pictures I crop the large picture (4800x4800) into multiple pictures with a resolution of 320x320. To identify each of those cropped pictures I calculated its center coordinate as follows (from upper left to right then to south):

dlat = (side * 360) / (2 * np.pi * r) #dlat for 32m in degrees
side = 32
r = 6371000 #avg earth radius

minx = float(minx)
miny = float(miny)
maxx = float(maxx)
maxy = float(maxy)

identifier = (minx, miny, maxx, maxy)

# Takes a 4800x4800 image tile and returns a list of 320x320 pixel images

tile = np.array(tile)

images = []
coords = []

N = 0
S = 4800
W = 0
E = 4800

# y coordinate is dlat/2 degrees, i.e. 16 meters, south of the maximum y coordinate.
y_coord = maxy - dlat/2

while N < S:

    W = 0

    x_coord = minx + (((side * 360) / (2 * np.pi * r * np.cos(np.deg2rad(y_coord))))/2) #16 m to the middle

    while W < E:
        # The first image is taken from the upper left corner, we then slide from left
        # to right and from top to bottom

        images.append(tile[N:N + 320, W:W + 320])
        coords.append((x_coord, y_coord))

        x_coord += (((side * 360) / (2 * np.pi * r * np.cos(np.deg2rad(y_coord)))))

        W = W + 320

    N = N + 320

    y_coord = y_coord - dlat

As I also use this approach to then convert the identified pixel polygons to real world coordinates I run into minor differences between the underlying map and the identified objects (10-20m difference) sometimes minor. So I ask myself what is wrong with my approach and where could I make it more efficient by using proper packages.

I actually did some research but it is rather confusing to me because I do not use GeoTIFF files or similar. I get a numpy array as output which I then rasterize with rasterio to retrieve the polygons in pixel coordinates. I convert these using the following function to real world coordinates (using the middle point of each picture as an identifier and fix point to calculate the other coordinates):

def centroid_coord(center_coord, distance_vector, size = 320):
'''
returns lat, lon of array centroid
:param center_coord: tuple of center coords, identifier of picture tile
:param distance_vector: tuple of distance vector, measured from center point in px
:param size: scaling factor for distance from center point in pixels (half the size of the segmented output (320)
:return: lat,lon array centroid
'''
dist_px_x, dist_px_y = distance_vector
x_center, y_center = center_coord
x_min = x_center - (((32 * 360) / (2 * np.pi * r * np.cos(np.deg2rad(y_center))))/2)
y_new = y_center + (side/size) * dist_px_y * (dlat/side)
x_new = x_min + (side/size) * ((size/2) + dist_px_x) * 360 * (1/(2 * np.pi * r * np.cos(np.deg2rad(y_new))))
return (x_new, y_new)

I think that this approach is not optimal and that there is a much more correct and efficient way using proper packages.

How can I do this?

Array: Binary array with a size (320,320) indicating where a certain object is or not. Using rasterio finding polygons of the objects in pixel coordinates.
Output: Convert polygons in pixel coordinates to world coordinates using the identifier of the input image (center point)
My problems:

Is the center point calculation (world-coordinatewise) correct like this?
Should I take the upper left corner point of the image to ease the process?
How can I convert the polygons with pixel coordinates correctly to world coordinate polygons?


texts: How to flag parcels as being within stop type dependent distance of transit stops?

texts:
first question here!

I have a parcel layer and a transit stop layer.  The transit stop layer has different bus stops, light rail stops, etc.

I want to add to my parcel layer attribute table fields that say if said parcel is w/in the buffer of each kind of stop.

the attribute table for my parcel layer would end up looking like

Parcelid  Muni Bart Bus OtherAttribute
1         1    0    1   189
2         1    1    0   130423

Where 1 is if the parcel is w/in the radius of the type of stop and 0 if it is not.

Is this sort of thing possible in ArcGIS, or do I have to write something...

Thanks!


texts: Populating attribute table with sequential numbers using ArcPy?

texts:
Using Python, how can we populate a field in a feature class automatically in increasing order?

The user must be able to provide the start number. For example, say he wants to populate values from 51 to 75 then he must be able to provide the start number that is 51.The user must provide the number after the script has been executed,for example if we have to populate the parcel_id in a fc named parcel then hows it possible,


texts: One to Many relation in ArcGIS for Desktop?

texts:
I'm working on ArcMap, and I have a one to many relation between a building column and the residents in the building.

The attribute table for the building has these columns: id,label,number.

I have an excel sheet containing these columns: number,resident

I'm trying to relate the number from building attribute table and the number in the excel sheet.

How to do that?


texts:
Would be glad if someone helps me do this. I'm trying to convert a raster .tif file to a .shp file for some analysis. I'm using a 1.5GB data file and have installed GDAL using Conda. But if I run the code it says the following error:

  RuntimeError: This is a BigTIFF file.  BigTIFF is not supported by
  this version of GDAL and libtiff.

What should I do to solve this error? Please explain as I'm totally new to Python environment!

import os
from osgeo import gdal, ogr
import sys
# this allows GDAL to throw Python Exceptions
gdal.UseExceptions()

#  get raster datasource

src_ds = gdal.Open("IND_ppp_2015_v2.tif")
print (src_ds)


texts: GRASS in QGIS - ' Cannot read raster map region' (running r.lake.xy)

texts:
I want to use the GRASS function r.lake.xy on a DEM that I have.

I load the DEM raster into QGIS, it is not georeferenced. Then I create a New Mapset. I set Projection Not defined. I either leave the Default GRASS Region, or set it to 90/-90/180/-180 (I've tried both). I'm told my Mapset (named DEM5) is successfuly created.
I open GRASS tools and select r.lake.xy.
The first strange thing is the dropdown for Name of input elevation raster map doesn't allow me to select either of the options ('DEM5' or 'Permanent'). So I type in DEM5.
I type in Seed point coordinates and Water level and name my output file.
I click run but I'm given the error message: Cannot read raster map region - pathname: Unable to open header file for raster map DEM5@DEM5.


texts: How to judge a point with Lat and Lon about whether it is in land or sea?

texts: I just want to judge a point with Lon/Lat about whether it is in land or sea in my front-end program. I want to konw how to do it with JavaScript. Any WebAPI can be used freely?


texts: Calculating Slope Between Pixels Through Time using ArcPy?

texts:
I have a 9 year weekly time series (~500 raster grids of equal cell size and extent).  I'm interested in obtaining the regression line slope between pixels (Imagine stacking all 500 grids on top of one another and running a linear regression between each individual pixel).  Essentially, this would measure the net change between pixels through my time series.

My main question is about the coefficients defined in the first answer to Making linear regression across multiple raster layers using ArcGIS Desktop?

Is the coefficient defined as 12 due to the temporal scale of the original question being at the monthly level and, if so, should my code's coefficient be set to 7 due to the weekly temporal scale between grids?

I have been referencing these Q&As that are similar:

Representing trend over time?
Making linear regression across multiple raster layers using ArcGIS Desktop?

I wrote this python code after reading the links referenced above:

from __future__ import division
import arcpy
import os
from arcpy import env
from arcpy.sa import*

# define workspace
arcpy.env.workspace=r"C:\Users\mbs7038\Documents\Remove_Gauge_Outliers\test\forSlope.gdb"

# enable overwriting
arcpy.env.overwriteOutput=True

# check spatial analyst extension
arcpy.CheckOutExtension('Spatial')

# define output paths
slopePath=r"C:\Users\mbs7038\Documents\Remove_Gauge_Outliers\SpatialTrends\SlopeTEST.gdb"

# list all rasters in the workspace
rasters=arcpy.ListRasters('*')
# sort rasters numerically
rasters.sort()

# get the number of rasters
n=len(rasters)
print n

# setup index
i=1

# define division
seed=(n*n*n)-n
print 'the global seed is {0}'.format(seed)

for raster in rasters:
    print i
    coef=(12*(i)-((6*n)+(6*1)))/seed

    print 'Raster {0} is {1}:'.format(i,raster)
    print 'the coef for raster {0} is {1}'.format(i,coef)

    # Multiple raster by coefficient
    if i==1:
        outSlope=(Raster(raster)*coef)
        i+=1  # same as saying i=i+1
    else:
        print 'adding {0} to outSlope'.format(raster)
        outSlope=outSlope+(Raster(raster)*coef)
        i+=1
    if i==6:
        break

# Save final slope grid
print 'saving final slope grid'
outSlope.save(slopePath + "/" + "floodChange_py")

print 'script is complete'

I created this code as a test which appeared to work on a 5 week subset of my data.
If what I have referenced above is correct then this code will work on a time series (of equal cell size and extent) of any length.  Although there were no bugs in the code, I have a feeling I may have made a mistake someplace.

On a side note I'm using PyScripter, hence the "from future import division" in the first line.  This disables the automatic floor function which rounds floating point numbers to the nearest integer away from zero.


texts:
I am trying to georeference a large map of iceland, and with every setting I try, the image streches out on the x axis. I have a feeling it has to do with coordnate projection, but I have tried everything I can think of (probably have tried at least 15 different crs). I have also tried every sampling method.
What else can I try to make it more of a square image like shown below?


texts: Filtering a WMS layer in QGIS Server doesn't accept "IN" operator

texts:
I am using OpenLayers to connects to QGIS Server, to serve a WMS layers.

I want to filter the WMS layer according to one attribute. I have followed the documentation   to do so. Everything is working fine if I use the '=' operator, like the following, for example:

FILTER=myLayerName: "myAttribute" = 1

But I need to have the criteria using the 'IN' operator, because I am filtering on multi values, like the following, for example (doesn't work)

FILTER=myLayerName: "myAttribute" IN (1, 2, 3)

Could you please advice me with the proper way to use the  'IN' operator with Integer values?


texts: Connecting geocoded point to an identify function

texts: I am using the esri javascript api to make a web map. I need to use an address locator to geocode a point and then depending on where that point lands information populates a side panel. I am assuming you have to pass the coordinates to an identify function but I have no idea how to do this. Basically, I want the person to enter an address, it zooms to the address and information about the area populates the side panel. The information in the side panel would come from attributes from various polygon layers on the map. Any sample codes would be appreciated. Especially how to take the geocoded result and populate a identify function with the coordinates. (Or however this is handled) Thanks


texts:
I would like to create a vector layer, add a point to it and then use it in qgis algorithms.

This is what I wrote :

#create temporary layer containing center point
    point=QgsFeature()
    point.setGeometry(QgsGeometry.fromPoint(QgsPoint(center[0],center[1])))
    pointlayer=QgsVectorLayer("Point?crs=epsg:2154", "temporary_points", "memory")
    pointlayer.startEditing()
    pointlayer.addFeatures([point])
    pointlayer.commitChanges()

    #intersection between point layer and polygon layers to get names codes
    fieldname_vector=processing.runalg('qgis:intersection', parcel, pointlayer, False, None)

But when the algorithm runs, it returns

    Unable to execute algorithm
Wrong parameter value: Point?crs=epsg:2154

Is there a way to use my layer without saving it, knowing that I'll erase it in the future ?


texts:
I installed GeoServer 2.12.1 on windows 8.1 x64 and I am trying to design a custom "login" form using AngularJS5 with Node.js, but when I make a request to geoserver it rejected by CORS, and get the error:

  Cross-Origin Request Blocked: The Same Origin Policy disallows reading
  the remote resource at http://localhost:8080/geoserver/rest/. (Reason:
  CORS header ‘Access-Control-Allow-Origin’ missing).

I read several questions here and follow the steps to enable CORS in geoserver by: downloading the matching "Jetty-Utility Servlets" jar file and putting it in the path:

C:\Program Files (x86)\GeoServer 2.12.1\webapps\geoserver\WEB-INF\lib

and also I uncomment the filters in the web.xml file and then restart geoserver .. but all that not worked with me.

I don't know what is wrong. As I said above I read several post about the problem but no luck.Like this question: Enabling CORS in GeoServer (jetty)?

The Solution:
I used in my request headers, the Authorization header which doesn't allowed by default by CORS, so I added it explicitly in the web.xml file:

<init-param>
<param-name>allowedHeaders</param-name>
<param-value>Content-Type,X-Requested-With,accept,Origin,Access-Control-Request-Method,Access-Control-Request-Headers,Authorization</param-value>
</init-param>


texts: Projection issues with QGIS 3.16 Swellendam Training data

texts:
I am following this QGIS 3.16 training, following it exactly as shown (https://docs.qgis.org/3.16/en/docs/training_manual/vector_analysis/basic_analysis.html), using their provided training data. And I am having issues with the projection. . I tried deleting the original data, then re-downloading it, and still I am getting this issue. So I don't know if there is something wrong with the CRS, or am I doing something wrong with the Buffer analysis (although it's straight forwards so don't see where I am going wrong there)?
Any ideas?
EDIT: I have tried saving as a shapefile and working through the same process, and I still get the same issue with a ? next to the layer. And subsequently, a mal buffer analysis.
I then tried uninstalling QGIS 3.16, and downloading 3.10, and I get the same issue.


texts: Apply Symbology to Layer not working outside of ModelBuilder?

texts:
I am running ArcGIS 10.2.2 Version 10.2.2.3552, using ArcMap.

I am creating a python script and one of the things i need to do is apply symbology from layer to a shapefile in my mxd that is maintained outside of the mxd.  It works fine in model builder but when i run it in python i get multiple errors.

Code:

# Import arcpy module
import arcpy

# Local variables:
Tioga = "Tioga\\Tioga"
Tioga_lyr = "C:\\Users\\jcarl\\Desktop\\Dynamic LH Map\\Symbology\\Tioga.lyr"

# Process: Apply Symbology From Layer
arcpy.ApplySymbologyFromLayer_management(Tioga, Tioga_lyr)

Runtime error  Traceback (most recent call last):   File "", line 10, in    File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 6004, in ApplySymbologyFromLayer     raise e ExecuteError: ERROR 000966: The input layer is not valid Failed to execute (ApplySymbologyFromLayer).