texts: Rasterio.mask returns error 'shapes are outside bounds of raster' - but they are definitely within raster bounds

texts:
I'm using rasterio.mask.mask to mask a raster with a shapefile.

import rasterio
import fiona

image_object = rasterio.open('/home/cate/image_name.img')

aoi_shapefile = fiona.open('/home/cate/AOI.shp', "r")
aoi = [polygon["geometry"] for polygon in aoi_shapefile]

image_aoi, image_aoi_transform = rasterio.mask.mask(image_object, aoi)

and I get UserWarning: shapes are outside bounds of raster. Are they in different coordinate reference systems?

As far as I can tell the shape is definitely within the raster bounds..

image_object.bounds
BoundingBox(left=117.4095061135422, bottom=33.73209332058161, right=120.49755473423146, top=35.64128279392083)

aoi_shapefile.bounds
(119.82677788035676, 34.076577156414714, 120.37042708436398, 34.575927115244426)

image_object.crs
CRS({'init': 'epsg:4326'})

aoi_shapefile.crs
{'init': 'epsg:4326'}

Here is the view in QGIS:

I'm using rasterio 1.0a10 - any ideas what is happening here?


texts: Need dojo array.map to wait for feature layer query to complete

texts:
I am executing an IdentifyTask and then using dojo array.map to create an array of Features to pass to an InfoWindow.  The stripped down version of my code follows.

        var deferred = identifyTask.execute(identifyParkParams);
        deferred.then(function (response) {
            featureArray = array.map(response, function (result) {
               if(test1){
                  testTemplate = new InfoTemplate(stuff1);
               } else {
                  featureLayer.selectFeatures(query, function(featureSet){
                    //processing code
                      testTemplate = new InfoTemplate(stuff2);
                      }
               }
            feature.setInfoTemplate(testTemplate);
            return feature;
                });
                return featureArray;
            }).then(function(featureArray){
                    app.map.infoWindow.setFeatures(featureArray);
            });

My issue is if test1 is false the array.map puts a null testTemplate in the feature instead of waiting until the selectFeatures is done processing.


texts: Failure creating output file: (OGR error: file location not a directory)

texts:
I am using the MMQGIS plugin to join sets of data on QGIS. I have been failing to do this for days on account of the output file location being rejected.
This error message appears:
Failure creating output file: Creation of data source failed (OGR error: C:\Users\cavca-admin\Documents\imd.shp is not a directory.)
It obviously suggests that I need to be choosing a directory as the output file location but I am unable to do this because it isn't valid. I have tried saving all relevant files in the one directory and choosing a range of file locations but nothing has worked.


texts: How to classify different surfaces having almost similar Reflectance responses?

texts:
I am using ENVI 5.3 to classify a Region of Interest clipped from Landsat 8 OLI. I have already performed calibration of my DNs by transforming it into Radiance and derived Reflectance image using FLAASH atmospheric correction module provided in ENVI. Afterward, I have re-scaled the output (Reflectance Values) to 1-0 scale prior to extract landcover classes from the scene.

Apparently, their is minute difference of Reflectance between builtup/impervious surfaces and exposed bedrocks. Using, Maximum Liklihood Classification doesn't seems to be working for the reason. It places both materials in a single class.

Manually editing the derived landcover classes like 'recoding' from vectors or masking such pixels prior to classification and afterward combine masked pixels to desired class is hectic task when we are dealing with large areas having similar spectral issues.

Kindly have a look on the below scene and spectral profiles of the both materials and suggest me how could I deal with these pixels having almost similar reflectance responses. What is the optimal way to classify such scene?


texts: Getting list of layer names using PyQGIS?

texts:
I need a list containing the names of all the layers in a QGIS session. I did the task as

layersNames = []
for i in self.iface.mapCanvas().layers():
   layersNames.append(str(i.name()))

but this has the problem that only the names for the visible layers are extracted. How can I get a list with the names of all (visible or not) layers using PyQGIS?


texts:
Last week I try unsuccesfully to install the plugin to open ecw files in QGIS 3.2.on Ubuntu 18.4.
Since then GDAL algorithms are not working for rasters. I have the exact same error for all of them:

Proces algorithm:
Algoritmo 'Orientación' comenzando…
Input parameters:
{ 'BAND' : 1, 'COMPUTE_EDGES' : False, 'INPUT' : '/home/basi/Escritorio/GEA Forestal/Proyectos/Biomasa/Mapas/MDT/MDT Gu.tiff', 'OPTIONS' : '', 'OUTPUT' : '/tmp/processing_32b112e632e44e699178c61428af21c2/bd7b41926c554ce5811f8c5bc6954822/OUTPUT.tif', 'TRIG_ANGLE' : True, 'ZERO_FLAT' : True, 'ZEVENBERGEN' : False }
GDAL command:
gdaldem aspect "/home/basi/Escritorio/GEA Forestal/Proyectos/Biomasa/Mapas/MDT/MDT Gu.tiff" /tmp/processing_32b112e632e44e699178c61428af21c2/bd7b41926c554ce5811f8c5bc6954822/OUTPUT.tif -of GTiff -b 1 -trigonometric -zero_for_flat
GDAL command output:
gdaldem: error while loading shared libraries: libNCSEcw.so.5.3.0: cannot open shared object file: No such file or directory
Execution completed in 0.05 seconds
Results:
{'OUTPUT': }
Charging result layers
Next layers were not generated correctly./tmp/processing_32b112e632e44e699178c61428af21c2/bd7b41926c554ce5811f8c5bc6954822/OUTPUT.tif

I tried to reinstall gdal and qgis, but same problem persists. Any idea?


texts: Does transforming to a new projection, then back, affect data accuracy?

texts:
I have a feature class (South Carolina counties, so a fairly large geographic area) in NAD83 SC State Plane. It needs to be transformed to a second projection (NAD83 UTM 17), then transformed back to the original. I'll be using Esri's Project tool to accomplish this.

Can this dual transformation cause a shift in the location of the polygons coordinates, and by how much -- centimeters, meters, kilometers?


texts:
I am trying to visualize the trajectory, but the online Tianditu map cannot be loaded.
I'm using QGIS 3.14.16-Pi, Python3.7.
def remoteAction(self):
    urlwithparams = 'https://t0.tianditu.gov.cn/img_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&' \
                    'LAYER=img&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}&' \
                    'tk=1d28355a0563c9409770728a2435940f'
    print(urlwithparams)
    tdtImgFull = 'type=xyz&url={0}&zmax=19&zmin=0&crs=EPSG3857'.format(
        str(QUrl.toPercentEncoding(urlwithparams, b":/?="), encoding="utf8"))
    print(tdtImgFull)
    rasterLyr = QgsRasterLayer(tdtImgFull, 'tianditu', 'wms')
    if rasterLyr.isValid():
        QgsProject.instance().addMapLayer(rasterLyr)
    else:
        print('invalid layer')

This is the output of my program：
https://t0.tianditu.gov.cn/img_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=img&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}&tk=1d28355a0563c9409770728a2435940f
type=xyz&url=https://t0.tianditu.gov.cn/img_w/wmts?SERVICE=WMTS%26REQUEST=GetTile%26VERSION=1.0.0%26LAYER=img%26STYLE=default%26TILEMATRIXSET=w%26FORMAT=tiles%26TILEMATRIX=%7Bz%7D%26TILEROW=%7By%7D%26TILECOL=%7Bx%7D%26tk=1d28355a0563c9409770728a2435940f&zmax=19&zmin=0&crs=EPSG3857
invalid layer


texts: Google fusion tables: Fill in lon/lat columns via geocoded address

texts:
I have a list of 100 companies that I want to geocode. I have all of their addresses but none of their lon/lat coordinates.

I have three columns related to location:

Address
Longitude
Latitude

Is it possible that when file->geocode is clicked (relating to the address column) that it fills in the longitude and latitude columns with the coordinates? I can then read these columns from there on in avoiding the need to geocode again.

Update

On further research it's seems I have to store the geocoded result in a database... not fusion tables as it doesn't seem possible with how I want it to work.


texts: Extracting multiple raster data to point shapefile in ArcGIS ModelBuilder?

texts:
I used raster iterator tool in ArcGIS ModelBuilder read several raster files in IMG format and extract values of corresponding raster files to a point Shapefile using "Extract multiple values" tool in spatial analyst. However, names of the raster files are too long and in extracted Shapefile fields they do not make any scene.

How can I extract part of the raster files name and use it as corresponding field name in the shapefile once the data are extracted?


texts: Difference between geometry and $geometry QGIS

texts:
This may be trivial question, but I'm confused, what is the difference between geometry and $geometry in QGIS? I didn't find any examples of using this function.

Particularly I'm interested in geometry, in what cases it is helpful?


texts: How to create Map from jpeg in Arcgis 10

texts:
I'm pretty new to Arcgis software and going to create map.
I have some points which should be shown on the map( there are temperatures of cites).
So, I need to create a map in arcgis. I guess!
What I have is a jpeg file of the area. Is there any way to do that?


texts:
I'm trying to run an NDVI calculation on a set of raster layers. The approach is to adapt this QGIS 2.x code to run in GQIS 3.x. A far as I'm aware it's the "legendInterface" that's changed, but I'm getting lost trying to suss out how it changed, and what to use as an alternative.
lddLrs = qgis.utils.iface.legendInterface().layers()
path = "/QGIS/Project/"
for lyr in lddLrs:
    entries = []
    #Define red band#
    ras1 = QgsRasterCalculatorEntry()
    ras1.ref = 'layer.name() +@121'
    ras1.raster = lyr
    ras1.bandNumber = 121
    entries.append(ras1)
    #Define nir band#
    ras2 = QgsRasterCalculatorEntry()
    ras2.ref = 'layer.name() +@192'
    ras2.raster = lyr
    ras2.bandNumber = 192
    entries.append( ras2 )
    calc = QgsRasterCalculator( '("layer.name() +@192" -  "layer.name()
   +@121") / ("layer.name() +@192" + "layer.name() +@121")', path + lyr.name()
   + "_NDVI.tif", 'GTiff', lyr.extent(), lyr.width(), lyr.height(),
   entries )
   calc.processCalculation()


texts:
I have a running Qgis Server 2.14 with a QGIS Web Client. As soon as I create a print layout in the composer and save the project in QGIS, it is no longer possible to load the map in the QGIS Web Client. It is stuck on loading map. This is the error message I get form the browser console.

Uncaught TypeError: Cannot read property 'width' of undefined
    at postLoading (WebgisInit.js:469)
    at Ext.tree.AsyncTreeNode.load (WebgisInit.js:246)
    at EXTUTIL.Event.fire (ext-all-debug-w-comments.js:518)
    at Ext.tree.AsyncTreeNode.fireEvent (ext-all-debug-w-comments.js:159)
    at Ext.tree.AsyncTreeNode.fireEvent (ext-all-debug-w-comments.js:54140)
    at Ext.tree.AsyncTreeNode.loadComplete (ext-all-debug-w-comments.js:55372)
    at Ext.tree.AsyncTreeNode.<anonymous> (ext-base-debug.js:969)
    at QGIS.WMSCapabilitiesLoader.processResponse (QGISExtensions.js:298)
    at QGIS.WMSCapabilitiesLoader.handleResponse (ext-all-debug-w-comments.js:56415)
    at Ext.data.Connection.handleResponse (ext-all-debug-w-comments.js:6528)

The code in WebgisInit.js:469 is this:

var mapWidth = composerTemplate.map.width / ptTomm;

The layout I set up in the composer is regular A4 300dpi


texts: Using ModelBuilder or Python parser of Field Calculator to check values between two datasets?

texts:
I have a *.dbf that I need to check values against for a shapefile.

I need to report whether the values match the dbf in the shapefile and, if not, what the value should have been (in two new fields)

Here's what I have:

DBF (correct values for all unique combinations of two fields reported in a third field), sort of like a look up table.
For each unique value within field A and B, there is a value in field C.
A shapefile to check the value in field C against the dbf based off fields A & B within that shapefile.

I know a very inefficient work-around (concatenating fields then manually joining them), but I'm trying to get more familiar with ModelBuilder and Python.

I prefer using ModelBuilder so if this can be done in ModelBuilder with maybe just a Python field calculate (simple if then statement) that would be ideal.

I have very little Python knowledge and limited ModelBuilder knowledge (basically I know how to drag and drop geoprocessing tools and set up parameters, but I struggle with the iterators and loops).


texts: Using update cursor to update feature class field from another table with incremental date/time query?

texts:
I have a feature class called AO15_16 where the field id needs to be updated.  The update will contain a unique id based on a query of another table where I am querying issueddate, issuedtime and `zipcode.  The ids are already in the the 2nd table based on a unique date, zip code and 1 hour increment per zip code.

The table below shows the hour increments from 6pm-7pm per zip code on the date of 7/9/2016

500 7/9/2016    6:42:00 PM  7:06:00 PM  0:24    GARDEN GROVE                92840   11  Variable    5.8 3707, 3709, 3736

501 7/9/2016    7:00:00 PM              92840       Variable    5.8

502 7/9/2016    6:42:00 PM  7:06:00 PM  0:24    GARDEN GROVE                92841   11  Variable    5.8 3707, 3709, 3736

503 7/9/2016    7:00:00 PM              92841       Variable    5.8

504 7/9/2016    6:42:00 PM  7:06:00 PM  0:24    GARDEN GROVE                92844   11  Variable    5.8 3707, 3709, 3736

505 7/9/2016    7:00:00 PM              92844       Variable    5.8

The table below shows unique rows for 7/10/16

506 7/10/2016   3:13:00 PM  5:56:00 PM  2:43    GARDEN GROVE                92840   25  SSW 8.1 3707, 3736, 3915, 3918

507 7/10/2016   4:00:00 PM              92840       S   9.2

508 7/10/2016   5:00:00 PM              92840       S   9.2

509 7/10/2016   6:00:00 PM              92840       S   8.1

510 7/10/2016   3:13:00 PM  5:56:00 PM  2:43    GARDEN GROVE                92841   25  SSW 8.1 3707, 3736, 3915, 3918

511 7/10/2016   4:00:00 PM              92841       S   9.2

512 7/10/2016   5:00:00 PM              92841       S   9.2

513 7/10/2016   6:00:00 PM              92841       S   8.1

The feature class contains an issueddate, issuedtime and zipcode which meets the conditions of the above table.

The goal is to update the field id in the feature class by reading the query from the table.  How can this be done with the query programatically changing so it takes into account the next greatest date and hour increment per zip code?

i.e. where issueddate = n + 1 AND issuedtime = n + 1 BETWEEN the hours that are in the rows?

What I have to use the table as a lookup table is below, how can this be changed to increment next greatest date and time in the where clause and update the 'id' field in the feature class?

whereClause = '''WHERE ISSUEDTIME BETWEEN x AND x (increment) and DATEADD(d,ROW_NUMBER() OVER (ORDER BY ID ASC), GETDATE())'''

 tableSet= set()
            with arcpy.da.SearchCursor(table, fields="Time Start; date; zipcode") as cursor:
                for row in cursor:
                    tableSet= .add(row[0])

            with arcpy.da.UpdateCursor(fc, "id", whereClause) as cursor:
                for row in cursor:
                    if row[0] in tableSet:
                        row[1] = row[0]
                    cursor.updateRow(row)


texts: I have an address locator that is updated weekly and I need to overwrite the existing geocode service and share the updated address locator as a new geocode service. Is there a pythonic way to automate the overwriting and re-publishing step, like in the screenshot below?


texts: Spatial join for two polygon files in QGIS?

texts:
I have two polygon files. One is MSOA (area defined by UK Census) area and the other is area defined by a project purpose. Most of the project area/polygons are within the MSOA areas e.g. five MSOA areas go within one project area. This happens in most cases. But there are few cases where MSOA areas are divided by project areas.

I want to join these two polygons. So that I can have a table saying MSOA area 1 covers project area 1,2,3 and 5.

Now I understand there will be at least two different methods to follow. I did this in Mapinfo where the join command works with two/three different options such as 'within' 'consists' etc. Unfortunately, I don't have Mapinfo atm and only have access to QGIS. The attribute tables do not have anything common. So it needs to be a spatial join.


texts: I am currently doing a QGIS project for my apprenticeship and I'm trying to understand the best way to get borehole log data into QGIS so I can visualise it and analyse the data.


texts: How to have "Preserve Layer and Table IDs" option enabled by default

texts:
A ability to preserve layer and table IDs between map service updates is a great feature added at 10.3.1. Unfortunatly it seems that it has to be manually applied to each dataframe when an .mxd is created.

Is there a way to have it enabled by default? Map publishers rarely remember to enable it and I occasionally end up having to edit my applications after updates anyway.

Edit
Found a relevant thread on geonet, it looks like you can't set this option through ArcObjects but you can set/change layerIDs on individual layers

Edit 2
Link to an idea I submitted at the ArcGIS ideas web site in case anybody else would find this useful.


texts: export more than one but not all attributes with qgis2web

texts:
Using the qgis2web plugin:
Is it possible to select more than one but not all attributes to be shown as popup?
I can select all attributes or only one attribute to be shown. but not for example two attributes.
Neither shift-click nor right-click seem to work. Any workaround (beside deleting all not to be shown attributes and then select all)?


texts: QGIS Server Accesing wms/wfs layers using user/pass

texts:
I have seen that there is an authentication system in QGIS which, after reading it two times, I don't understand very well:

https://docs.qgis.org/2.14/es/docs/user_manual/auth_system/auth_overview.html

I am trying to publish in QGIS Server some postgis layers but they have to be accesible only using a user/pass. I think that what I am reading in that article has nothing to do with what I am trying to achieve.

Is this possible using QGIS Server/Client architecture?


texts: Extracting point data from EPSG 4326 image

texts:
I am new in the field.

I am trying to extract remote sensing data from selected points. I am using a remote sensing image with CRS of EPSG:4326. I defined the points which I will extract data from by entering their longitude and latitude as decimal degrees. (I was given this information both in decimal degrees and in degrees - minutes - seconds. I checked them using Google Maps to see if they correspond to my locations of interest.) I am using QGIS for extracting the point data from remote sensing image.

The question I have is; since the EPSG 4326 is not a projection, would the point data I get from the image be accurate this way? Or should I project the image first before extracting the data?

I am only interested in the remote sensing data I will extract from the point locations.


texts: Copy Parallel editor tool - python script / modelbuilder

texts: Is it possible to use Copy Parallel tool in a python script or in modelbuilder for ArcGIS 10.2? Does anyone have a workaround or script for a tool that essential does the same thing as Copy Parallel. I have a modelbuilder model that needs to copy parallel lines from a street centerline shapefile, buffering will not work in this case.


texts: How to move the QGIS Layerpanel/Dockwidget programmatically

texts:
i'm working on a QGIS-plugin and want to move the default LayerPanel/Dockwidget to somewhere else programmatically. More specific: i want to move it into a tab of my tabwidget.
Is this possible? How can i get access to that specic dockwidget out of python?


texts: Why the modifications I do in QGIS doesn't appear when I save an OSM file I created?

texts:
Through QGIS I loaded an OSM file I created and I'm trying to visualize what happens with the file in terms of codes after some modifications like adding points, lines or polygons. But the problem is that when I save the file with the OpenStreetMaps plugin the modifications don't appear as I was expecting...

Here you can see the file I created and if you want, feel free to do the same test adding points or whatever, saving and then see if something changes in the new saved file.

<?xml version="1.0" encoding="UTF-8"?>
<osm version="0.6" generator="OpenStreetMap server">
<bounds minlat="52.2615973629" minlon="10.5262132592" maxlat="52.2653605629" maxlon="10.5313793802"/>
<node id="100" lat="52.264859" lon="10.528042" version="1" user="ribe_lu" visible="true" timestamp="2012-09-10T13:47:15Z"/>
<node id="101" lat="52.264859" lon="10.529796" version="1" user="ribe_lu" visible="true" timestamp="2012-09-10T13:47:15Z"/>
<node id="102" lat="52.262619" lon="10.528042" version="1" user="ribe_lu" visible="true" timestamp="2012-09-10T13:47:15Z"/>
<node id="103" lat="52.262619" lon="10.529796" version="1" user="ribe_lu" visible="true" timestamp="2012-09-10T13:47:15Z">
    <tag k="made_by" v="Lucas"/>
</node>
<node id="104" lat="52.264186" lon="10.529078" version="1" user="ribe_lu" visible="true" timestamp="2012-09-10T14:20:17Z"/>
<node id="105" lat="52.26326" lon="10.529078" version="1" user="ribe_lu" visible="true" timestamp="2012-09-10T14:20:17Z"/>
<node id="106" lat="52.26326" lon="10.528801" version="1" user="ribe_lu" visible="true" timestamp="2012-09-10T14:20:17Z">
    <tag k="made_by" v="Lucas"/>
</node>
<way id="25" visible="true" timestamp="2011-01-09T14:05:48Z" version="1" user="ribe_lu">
    <nd ref="100"/>
    <nd ref="101"/>
    <nd ref="103"/>
    <nd ref="102"/>
    <nd ref="100"/>
    <tag k="area" v="yes"/>
    <tag k="Schloss" v="Arkaden"/>
    <tag k="Postal_Code" v="38100"/>
</way>
<way id="26" visible="true" timestamp="2011-01-09T14:05:48Z" version="1" user="ribe_lu">
    <nd ref="104"/>
    <nd ref="105"/>
    <nd ref="106"/>
    <nd ref="104"/>
    <tag k="Test_by" v="Lucas"/>
</way>
<relation id="75" visible="true" timestamp="2012-09-10T14:50:10Z" version="1" user="ribe_lu">
    <member type="way" ref="25" role="outer"/>
    <member type="way" ref="26" role="inner"/>
    <tag k="type" v="multipolygon"/>
</relation>
</osm>


texts: Remove duplicate strings from a field with Field Calculator in ArcGIS

texts:
I tried using the same codeblock in this post:
Remove duplicates from a field, but it is just copying the contents and not removing duplicates. I think my data is formatted the same. Anyone see what I'm missing?

Codeblock

def delDups(n):
return ', '.join(set(n.split(', ')))

Expression

delDups( !Names3! )


texts: How to create a raster tiff (dem / heightmap) from a mesh of xyz points?

texts:
I have to create a terrain heightmap(raster tif preferably) from mesh of xyz ungridded points(in geographical coordinates), I tried taking min and max of bbox of mesh and making width & height by subtracting max-min and then querying for height on all the vertices in the mesh. But it works only if they are in geographical cordinate system n not lat-long as then it can be really tiny. Any ideas how to do it?

My code snippet to create Raster :

double startX = mesh.verticesBBox.bbmin.X;
double startY = mesh.verticesBBox.bbmin.Y;
double bottomX = mesh.verticesBBox.bbmax.X;
double bottomY = mesh.verticesBBox.bbmax.Y;

int width = fabs(bottomX - startX);
int height = fabs(bottomY - startY);

const int channel = 1;  //Dem
const char *pszFormat = "GTiff";
GDALDriver *poDriver;
char **papszMetadata;
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
GDALDataset *poDstDs;
int count = 0, a =0 , b=0;
poDstDs = poDriver->Create(tfile.toString().c_str(), width, height, channel, GDT_Float32, NULL);
GDALRasterBand* band = poDstDs->GetRasterBand(1);
std::vector<unsigned char> data(width*height*GDALGetDataTypeSize(GDT_Float32) / 8, 0);
int bufferSizeX = width*GDALGetDataTypeSize(GDT_Float32) / 8;
int bufferSizeY = height*GDALGetDataTypeSize(GDT_Float32) / 8;

for (///////////vertex array iterator)
{
    double h;
    dem.GetHeight(Vector2<REAL>(ver->X, ver->Y), h);
    data[count] = unsigned char(h);
    int locationX = fabs(ver->X - startX);
    int locationY = fabs(ver->Y - startY);
    if ((width > locationX) && (height > locationY))
    {
            CPLErr e = band->RasterIO(GF_Write, locationX, locationY, 1, 1, &data[count], 1, 1, GDT_Float32, 0, 0);
    }
        count++;
}


texts: Giving generated random raster same extent and coordinate system of input raster in ModelBuilder?

texts: In ArcGIS's (10.1, advanced license) ModelBuilder, I want to create a random raster (let's call it B) with values between 0.90 and 1.0, and I want to multiply another raster (let's call it A) by B. A is fed into the model by the user. My goal: I want B to have the same extent and coordinate system as A. What I think to do to accomplish the first goal is to drag and drop the CREATE RANDOM RASTER tool into the ModelBuilder and then using the tool's drop-down menu "output extent" to select raster A. I do not know how to proceed about the second goal, namely selecting for raster B the same coordinate system as A. I have added a picture portraying an illustrative example of the model.


texts: Unknown coordinate system need in Lat/Long degree format

texts: I received CAD data from PDs and the data they gave is already geocoded (which is great) and it reads well when I add it to my map on ARC map. However, the format is as follows (x = 632086.7, y = 1043211.35) I need them in the following format (x = 33.867022, y = -112.1409750) so that I can use it in R. My first question is what is the original coordinate system format? I've tried identifying it but I haven't found something that is helpful. Second question is, how can I transform the first coordinates into the second format?


texts: ArcObjects ITextSymbol.GetTextSize in Python: Why are the out params equal to 0?

texts:
This is related to the question I asked earlier about the syntax for this method under Python and comtypes.  Now that I know the syntax and am trying to use it, the returned out params are equal to 0.  Why is this?  Here is my code:

#make TextElement
pTextElement = NewObj(esriCarto.TextElement, esriCarto.ITextElement)
pTextFont = NewObj(stdole.StdFont, stdole.IFontDisp)
pTextFont.Name = 'Arial'
pTextFont.Size = 34.0
pTextSymbol = NewObj(esriDisplay.TextSymbol, esriDisplay.IFormattedTextSymbol)
pTextSymbol.Font = pTextFont
pTextSymbol.HorizontalAlignment = esriDisplay.esriTHACenter
pTextSymbol.VerticalAlignment = esriDisplay.esriTVACenter
txt = 'T>'   #((( hard coded for testing )))
pScrDisp = mxDoc.ActiveView.ScreenDisplay
pDisplay = CType(pScrDisp, esriDisplay.IDisplay)
#Get dimensions of text string
pDisplay.StartDrawing(0, -1)
dblX, dblY = pTextSymbol.GetTextSize(pDisplay.hDC, pDisplay.DisplayTransformation, txt)
pDisplay.FinishDrawing

The out parameters dblX and dblY both equal 0.  The hDC and DisplayTransformation inputs are valid.  I checked.  So, there is something I'm missing.  Any insights?


texts: Searching on a remote GeoNetwork catalogue

texts:
On a GeoNetwork 3.0 node, I would like to search directly on another remote GeoNetwork node rather then harvest its data.

My goal is to incorporate search results from the remote node and the local node into one single response.

I can search on a local catalogue using the q-search service from the API.

In the docs from version 2.10 I found the hint, that there is an option remote, which should be used for searching on remote GeoNetwork catalogues when using the xml.search service.

I wonder, if this can be applied for the q-search service in the 3.x version. Is there any example how to achieve this or has anyone done something similar before?


texts: Calculate PV rooftop potentials in Berlin QGIS 3.16. (OS Windows 10)

texts:
I want to analyse PV potentials of roofing in QGIS (v. 3.16).
I am experienced in vector data analysis but have barely worked with raster data sets. I therefore have no experience and am asking for assistance with the task.
I have access to the digital surface model of Berlin (as WMS layer) through the following link:
https://fbinter.stadt-berlin.de/fb/wms/senstadt/k_dom
I have read previous posts about similar analysis like:
Calculating total PV potential of roof using QGIS
and am wondering if someone could guide me through a similar analysis.
The SAGA toolbox that I can open in QGIS does not offer the Potential Incoming Solar Radiation algorithm.
I have further installed a SAGA standalone (version 7.8.1) that offers the toolkit.
Is there another way to calculate Potential Solar Radiation?
Side note: I have installed the following programs on my computer (OS Windows 10):

QGIS (3.16.1.) standalone
SAGA (3.8.1.) standalone
WinGRASS (7.8.4) standalone
Python (3.8.6)


texts: Problem accessing 10.0 File Geodatabase from ArcGIS Desktop 10.1

texts:
I'm on 10.0 SP4 (build 4000) and two co-workers are on 10.1 original release (build 3035 I believe). The DB was provided by a client and saved on a network drive. Since I have no issue accessing the data I assumed it was created in 10.0 or earlier. One of the co-workers is trying to access the data remotely, the other from within the office network. They both got the same "Failed to connect to database" error message. A third co-worker on 10.0 gave it a try in the office and it worked for him.

This looks to me a typical backward compatibility scenario, plus there are other DBs created in 10.0 and they had no issue accessing from 10.1. I'm puzzled why this is happening? Shapefiles were created so everyone can reference the data, but I'm really curious in finding out what could be the cause of the issue.


texts: Should a GDAL Plugin be the same version as GDAL aready built into QGIS

texts:
Should a GDAL Plugin be the same version as the GDAL that's aready built into the QGIS version installed?

Or,can another different version of GDAL be built and it installed as a Plugin for QGIS, even though the version of GDAL is different?  Is it reading QGIS's GDAL installation or a standalone GDAL installation for the plugin?

I'm going to try and create a Ubuntu QGIS Lisboa plugin with File GDB support.  And, I'd like to do the same for Windows 7 QGIS Lisboa. It's going to be installed using a standard osgeo4w (rather than advanced) install, but I'd like to build a plugin for it.  The reason I'd like to use standard install is that I'd like to have all the standard Python, PostGIS and Spatialite packages installed as well, and not sure how to do that with "Advanced", and still get GDB support.

If I can't get QGIS to use a new GDAL plugin featuring GDB support, the I'd like to atleast be able to add another standalone version of GDAL, and use the command line ogr2ogr to convert File GDB to shapefiles and dbfs (for the additional domain tables). Thanks in advance for your advice.


texts: Importing file geodatabase (GDB) to table of contents (TOC) of ArcMap?

texts: I have a geodatabase with several feature dataset and feature classes associated with it respectively. I would like to create a .mxd where the Table of contents consist of group layer and each group layer has several shapefile. The group layer name should be same as Feature dataset name.The feature classes in each feature dataset comes shapefiles with in that particular group layer as shown in figure.


texts:
I am trying to import  CSV file and get no reported errors but there are no points visible. If I click on "zoom to extent" nothing happens (even if I move away from the location of the points before doing this.  I have shut off all raster layers but they were not hiding the points.  I have tried to find an answer here but none seem to fit my problem. I have followed the steps in this tutorial: http://www.qgistutorials.com/en/docs/importing_spreadsheets_csv.html

Here is an image of the CSV file:

here is the import dialog box:

and here is the QGIS attribute table after import:

I am saving in WGS84.  I tried the local projection (NAD83 Zone 19N) but that made no difference.

I am using QGIS 2.12.1 on Windows 10.
Any suggestions on what I am doing wrong?


texts: Openlayers dynamic viewport width and height assignment?

texts:
I display an openlayers map in a html page by placing it in a  tag and assigning a css class to it so as to set its dimensions (map size and map bounds).

i then attempt to dynamically create a style element (my custom division class), append it to the document and change the map class to this custom style.

2 unintended things happen:

just appending the style to the document changes the class of the map division (bad DOM on my part, surely);
the map division still uses "old" bounds (does not center correctly on zoom);

Can anyone recommend a better approach to this situation?

Code and further details on this stackoverflow post


texts: Is it possible to restore lost PRJ-Files?

texts:
Currently I am working on a project with QGIS (2.14.1). Therefore I need to add old (but important) shapes that were drawn by a colleague.

Somehow most of the .prj-files are missing - I don't know why, maybe they were deleted by accident or not created at all. Those shapes with mising prj meta-data are not displayed correctly, depending on the CRS I choose for them they eather appear on a random place on earth or are not displayed at all.

I don't have a clue what I can do, is it even possible to fix that somehow? Maybe via reprojection or something?

I don't know the original projection btw.

The project area covers parts of Southwestern-Germany, County of Rhineland Palatinate, close to Mainz. Based on that I created prj files via spacialreference.org for potential CRS (eg. EPSG 31466, 31467, 3044, 31493, 3857, 32632). Unfortunately without success.

Layer extent:  xMin,yMin 3430186.73;5517809.87 : xMax,yMax
   3443672.61;5529390.81


texts:
In print composer I have one map frame with an overview of the main map frame.
It used to work- updating as I changed the scale of the main map.  But now when I open up qgis and choose the same map layout in print composer the overview frame does not show.  Other map layers draw in the same map frame just fine but the overview frame is missing.
I am using the most recent version of QGIS 3.16.1
If I add a new map frame and add an overview to the new map it works.  But the original map I was using for this purpose does not.
Thanks for any insight.


texts: I am trying to analyze the service area of bus stops in that yellow area. So basically those bus stops contain the bus stop id and their latitudes & longitudes. But then there's those bus stops that I circled in red that doesn't produce a polygon like the rest of the other bus stops. Is there any proper way on placing the points to properly produce service areas? Put in mind that i used the service area analysis in ORS Tools plugin


texts: Finding polygon from two arcs using CLI tools or some library

texts:
I'm trying to convert from textual representation into a set of arcs and lines.

The 'input' I have is:

The area is within lines joining two semicircles with a radius of 0.67 NM centred on:
641139N 0214540W and
641244N 0214130W.

What I need is to be able to output similar to:

ARC CenterPoint, Radius, Start, End
Line End, Start2
Arc CenterPoint2, Radius, Start2, End2
Line End2, Start

So, what I am looking for are the Start, End, Start2 and End2 points.

I have no problem writing small python programs etc to do this, but I'm looking for some libraries that can help me with the GIS part.

For clarification, the shape I'm trying to describe can be seen on page 'BIRK AD 2.24.11.1 - 1' in http://aip.samgongustofa.is/media/PDF/AD_2_BIRK.pdf


texts: How to fix label clipping in Mapbox Studio

texts: I have labels that are being cut off when rendered in Mapbox Studio. As I understand it there is a 256x256 tile that gets rendered and if the label text "bleeds" into another tile it may get cut off. Research led me to a 'buffer-size:' thing but that isn't making any difference whether I set it to 0, 512, or 60000.


texts: Deleting large numbers of fields in ArcGIS 10

texts:
It there a good strategy for deleting a few thousand fields from the attribute table of a shapefile in ArcGIS using either ArcMap or ArcCatalog?

I know how to manually right-click and choose 'delete field;; I'm wondering if there is any other strategy. If I could export the shapefile with columns selected, that would also be fine, but I don't think that can be done.


texts:
Is it possible to group rows together that have geometries that overlap?  If so, is it also possible to make that a configurable threshold.  Let's say I have a table representing data that looks like:

Let's say I'd like to get the number of geometries in each overlapping cluster together with the bounding box of the cluster.  I'd like to do something like:

SELECT
    ST_Extent(geom) as bbox,
    COUNT(*) as count
FROM table
GROUP BY ???

The result that I would hope for is:

|bbox     |count|
|---------|-----|
|bbox(B,C)|  2  |
|bbox(D,E)|  2  |
|bbox(A)  |  1  |

As a bonus question, I'd also like to be able to parameterize this over an overlap threshold like Jaccard index.  Meaning I'd like to get the same sort of grouping if the geometries overlap sufficiently.  For example, I'd like to get the same grouping as previously mentioned but split up D and E if the don't overlap by the supplied threshold:

|bbox     |count|
|---------|-----|
|bbox(B,C)|  2  |
|bbox(D)  |  1  |
|bbox(E)  |  1  |
|bbox(A)  |  1  |


texts: Is it possible to get all projections form OGR?

texts:
I'm working on a project where the user can give my program GIS data which may or may not have the projection it is using.

Is it possible to get a list of all the projections from OGR or Proj4?

Edit
The program I am working will take a GIS file from the user as an input and will show all the points. But sometimes the files do not say which projection they are using so I don't know what to do with the data. Instead of just throwing an error message I am trying to show a tree view of all possible projection the user can choose from. Kind of like ArcGIS, where you can change which projection is being used. Also I'm using c++ for this project.


texts: Why difference between Google Earth and QGIS distance measurement?

texts:
I am quite new to QGIS.

I wanted to test if Google Earth and QGIS distance measurement tools are actually getting the same results. Turns out, they don't -at least not in my case.

I'm pretty sure the project CRS (Stockholm, Sweden) is the right one - and I guess Google is always right (!)?

The CRS is SWEREF99 18 00, EPSG::3011 (OTF).

Method:

DRAWING LINE: new layer -> draw line
MEASURING in QGIS: attribute table -> field calc -> $Length/1000*
EXPORT TO GOOGLE EARTH: save layer as.. -> KML-file (aware of correct CRS)
MEASURING in GOOGLE EARTH: import KML-file -> measure manually w/ GEarth measurement tool

(please look at the picture to see how different they measure 3 different lines, *= I've seen this in several tutorials)

So, my question is either;

1) Which measurement should I trust - QGIS or Google EArth?

and/or

2) How do I get QGIS to measure correctly?


texts: qgis gdaltransform versus OSGeo.OSR.CoordinateTransformation

texts:
I have implemented a simple coordinate transformation test in c# using GDAL 2.4.1 downloaded from https://www.nuget.org/packages/GDAL/. When comparing the result with gdaltransform (QGIS 3.4) and EPSG.io I do encounter a small difference in result.

c# test:

SpatialReference spatialSrc = new SpatialReference(string.Empty);

spatialSrc.ImportFromEPSG(32028);

SpatialReference spatialDest = new SpatialReference(string.Empty);

spatialDest.ImportFromEPSG(4326);

CoordinateTransformation transformator = new CoordinateTransformation(spatialSrc, spatialDest);

double[] result = new double[3];

transformator.TransformPoint(result, 2397350, 386280, 0);

result[0]=-76.3056550480934

result[1]=41.2177849636383

QGIS gdaltransform test:

C:> gdaltransform -s_srs EPSG:32028 -t_srs EPSG:4326

2397350 386280

-76.3053143354412 41.2178680673729 0

Result from EPSG.io (identical to the QGIS result):
-76.3053142 41.217868

I seek any explanation of the difference between my c# test and the QGIS/EPSG.io test results.

Am I missing something essential?


texts: Can GDAL ReprojectImage run multicore/multithread from python?

texts:
I am using GDAL in python to reproject an image. Is there a way to use all CPU cores in this process. I have tried the following with no success:
gdal.ReprojectImage(vrt, output, input_proj, ref_proj, gdalconst.GRA_Average, options=['NUM_THREADS=ALL_CPUS'])

I also tried including the NUM_THREADS=ALL_CPUS option to the output dataset:
output = driver.Create(fn, ref_x, ref_y, 1, gdalconst.GDT_Float32, options=['NUM_THREADS=ALL_CPUS'])


texts: Attribute table for shapefile showing nonsense symbols on my PC

texts: I am working with a shapefile which has several hundred rows and 10 or so columns. Usually, the cells contain data, but at the moment they are just showing nonsense symbols (attached). I have tried using the same shapefile on a colleague's computer, and it worked fine. I have also reinstalled QGIS. What is wrong?


texts: Error while importing networkx module in ArcMap10.6

texts:
I defined a custom toolbox in ArcMap 10.6, the script ran perfectly when I first use  it, but ran in error when I use the same script with same parameters again. After I restart ArcMap 10.6, I can use it perfectly again.

When I use the same script in ArcGIS Pro,it run perfectly. I want to know if it is a bug? If so,  is it a ArcMap's bug or a networkx package's bug or others? And what can I do?

Here is the error message:

File "******\caculate_network.py", line 7, in <module>
    import networkx as nx
  File "C:\Python27\ArcGIS10.6\lib\site-packages\networkx\__init__.py", line 128, in <module>
    import networkx.drawing
  File "C:\Python27\ArcGIS10.6\lib\site-packages\networkx\drawing\__init__.py", line 6, in <module>
    from . import nx_pydot
  File "C:\Python27\ArcGIS10.6\lib\site-packages\networkx\drawing\nx_pydot.py", line 27, in <module>
    from pkg_resources import parse_version
  File "C:\Python27\ArcGIS10.6\lib\site-packages\pkg_resources\__init__.py", line 48, in <module>
    from pkg_resources.extern import six
  File "C:\Python27\ArcGIS10.6\lib\site-packages\pkg_resources\extern\__init__.py", line 45, in load_module
    mod = sys.modules[extant]
AttributeError: 'NoneType' object has no attribute 'modules'


texts: No outputs created using processing-TauDEM and QGIS 3

texts:
Has anyone run across this error when trying the New processing-TauDEM plugin in QGIS 3?

  The following layers were not correctly generated:

  C:/Users/david.laverty/AppData/Local/Temp/processing_3e96b8dc52a84af0a986623fa57e31f9/4262183f0f124676a16440b55b006d67/PIT_FILLED.tif

  You can check the 'Log Messages Panel' in QGIS main window to find
  more information about the execution of the algorithm.

My Log messages look like this:

  2018-10-25T11:31:04     INFO    C:/MicrosoftMPI/Bin\mpiexec -n 2
  C:/PROGRA~1/TauDEM/TauDEM5Exe/pitremove -z
  C:/Users/david.laverty/Downloads/Ardstraw_02_02_2012/Ardstraw_02_02_2012/DTM_1m/V0900583.asc
  -fel C:/Users/david.laverty/AppData/Local/Temp/processing_3e96b8dc52a84af0a986623fa57e31f9/4262183f0f124676a16440b55b006d67/PIT_FILLED.tif

I have tried this as a TEMP layer save and as a save layer option but still I get the same result. There is nothing in the Log Messages regarding TauDEM and in the TEMP folder where the output is supposed to be saved there is nothing.

My TauDEM is working correctly as I am able to run the commands from my QGIS 2.18 install and also from the command line.

Has anyone else had this happen to them using the new TauDEM provider plugin?


texts: Reclassify ESRI raster using variables from a SQL Server table (connecting with pyodbc)

texts:
I am trying to reclass a raster file (representing vegetation and land cover types) using values stored in a SQL Server Express table. I have a set of ~400 species records in the table, with each row representing a different species, and a column for each vegetation type. Each vegetation type is coded as either "1" if it is suitable habitat for that species, or "0" if it's unsuitable. The vegetation raster will then be reclassified to two values, "1" or "0", based on the values associated with the species record (this will be slightly different for each species).

In PythonWin, I'm using pyodbc to connect to the SQL Server Express database table, then executing a select query statement to collect the values for a species record (row) into the pyodbc cursor. Then I want to assign each column value to the output raster value in the remap statement (see code attached). Unfortunately, I keep getting the following error:

TypeError: list indices must be integers, not tuple

The cursor has all of the values from the SQL Server Express table based on the select query, and I can access the data through the cursor object's row parameter. So the row.BIOME_X variable is definitely storing the 1 or 0 value. Or am I just doing something dumb with the Python syntax. Any ideas?

Thanks!

# Load Python libraries
import pyodbc
import arcpy
from arcpy import env
from arcpy.sa import *
import os
arcpy.env.overwriteOutput = 1

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# set variables
modelList = arcpy.ListFiles() # build list of species model names for loop
biome_Cur = ("xxxxxxx/xxxx/xxxxxx/1_Input.gdb/biome_current") # the original raster which will be reclassed

# connect to SQL Server Express database
cnxn = pyodbc.connect('DRIVER={SQL Server Native Client 10.0};SERVER=DESKTOP\SQLEXPRESS;DATABASE=Species;UID=sa;PWD=XXXXXX')
cursor = cnxn.cursor()

# main processing loop
for model in modelList:

    # reclassifies biome raster based on suitability code from SQL Server Species database
    # select query
    cursor.execute("""
                SELECT [BIOME_1],[BIOME_6],[BIOME_7],[BIOME_8],[BIOME_9],[BIOME_10],[BIOME_14],[BIOME_19],[BIOME_20],
                [BIOME_21],[BIOME_22],[BIOME_23],[BIOME_24],[BIOME_25],[BIOME_27],[BIOME_29],[BIOME_30],[BIOME_31],
                [BIOME_32],[BIOME_35],[BIOME_36],[BIOME_37],[BIOME_38],[BIOME_39],[BIOME_40],[BIOME_41],[BIOME_42],
                [BIOME_43],[BIOME_44],[BIOME_45],[BIOME_46],[BIOME_47],[BIOME_48],[BIOME_50],[BIOME_100],[BIOME_200]
                FROM Species.dbo.BiomesPerSpp_Rehfeldt
                WHERE ID = ?""", (model))

    # assign remap variable for reclassification
    remap_cur = arcpy.sa.RemapValue([7, row.BIOME_7],[8, row.BIOME_8],[9, row.BIOME_9],[14, row.BIOME_14],[20, row.BIOME_20],
    [21, row.BIOME_21], [22, row.BIOME_22], [23, row.BIOME_23], [25, row.BIOME_25], [30, row.BIOME_30],
    [31,row.BIOME_31],[32, row.BIOME_32], [36, row.BIOME_36], [38, row.BIOME_38], [41, row.BIOME_41],
    [42, row.BIOME_42], [43, row.BIOME_43],[44, row.BIOME_44], [45, row.BIOME_45], [46, row.BIOME_46],
    [47, row.BIOME_47], [50, row.BIOME_50], [100, row.BIOME_100],[200, row.BIOME_200])

    biomeReClass_cur = arcpy.sa.Reclassify(biome_Cur, "Value", remap_cur, "NODATA")


texts: Is there a standard for oil pipeline symbology?

texts:
To develop a map for specific public is recommended using some special symbology. I am in charge to develop a map for an audience in oil pipeline business, I would like to produce my map using a standard symbology, but so far I have not found any.

I just read some post here that led me to the INSPIRE documentation, but it does not fulfil my requirements. Does anyone have information about this topic?


texts: Multi Ring Buffer "Outside polygons only" option grayed out - ArcGIS

texts: I'm trying to create a multi-ring buffer around a shapefile of ~400 points in ArcMap 10.6.1 using the Multiple Ring Buffer tool. I need the buffers to be "doughnut style" not "swiss cheese" style, where the larger rings may overlap each other or with smaller rings but remain separate annulus-shaped features. The "Outside polygons only" option for this tool becomes grayed out as soon as I enter the points shapefile I want to use for buffering. Is there a way around this where the option is not grayed out? I would like to be able to use this simple checkbox. Because the larger rings overlap with the smaller rings, I can't use Erase to create the rings either.


texts: IOError: Decoder Zip Not Found while using mapproxy with mapserver in windows 7 64bit machine

texts:
I am new to GIS technology.  I have configured a map server OSGEO4W.  I am able to render maps using WMS and openlayers. But the rendering is very slow.

For faster rendering I want to use Mapproxy.  I have installed Mapproxy with Python virtual environment.  Mapproxy service is running fine.  I am checking it with the following url "http://"mylocalhost":8080/"
Then, it is showing:

"Welcome to MapProxy 1.5.0"

I have changed the mapproxy.yaml as per the demo application file so that it shows its "demo" application. When I am browsing demo application and rendering map using EPSG:4326 or EPSG:900913 (even when I am using it in my application with openlayers and map file, then also it is showing the the following):

It is showing a blank image.
It is showing the following error on command window:

  File
  "c:\mapproxy_venv\lib\site-packages\pil-1.1.7-pu2.7-win-amd64.egg\Image.py"
  line 385, in _getdecoder raise IOError("decoder %s not available" %
  decoder_names) IOError: decoder zip not available

It means it does not receive "decoder zip".

Please help me to fix this problem.

Thanks


texts:
I've a created a web service using the ArcGIS API for Javascript that displays isochrones (service area polygons) generated from a network dataset built in ESRI's network analyst.

My issue is that the polygons are very jagged and I'd like to smooth them.  Typically, I accomplish this in the polygon generation tab of the service area properties in ArcMap (trim, generalized, detailed).

However, this seems to have no effect after publishing the dataset as a service.  I've been able to alter the polygons by adjusting the Default Output Geometry Precision both in the Service Editor and in the Javascript, but this alone is not enough.

Are there other options either in Javascript or in the publishing process for smoothing polygons that I've missed?


texts: QGIS rounds up the coordinates to six decimals when saving layer changes

texts:
I have the problem that when I create a new point (in a point layer) or a new line (in a line layer of course) and I save the changes to the layer (so NOT the project save) then QGIS rounds up the coordinates of the new element to six decimals. I can see that in the "identify results" window by clicking on "(Derived"). So as a result, I don't get the required position of the new element. The layers are DB tables that I received and I just open them in QGIS. I don't know the source of the layers. The CRS of the project is the same as the CRS of the layer. The thing is, if I create a new shapefile-layer myself (point for example), the new element remains in the wanted location even after saving layer changes. So it has something to do with the layer I received. But I don't know what that is. I am trying for 3 days to find a solution and can't seem to solve the problem. Does anyone have a hint for me?

Bevore saving the layer changes:

After saving the layer changes:


texts: GPS waypoints are giving me different coordinates once transformed into proper coordinate system?

texts:
I recently went and collected some waypoints with a GPS that is in NAD27 UTM16N. After uploading the GPX file into QGIS, I saved it as an ESRI shapefile, and changed the CRS from WGS84 to NAD27 UTM16N. After adding it to a new file, the coordinates appeared to be in the correct place, but after comparing them to the points in the actual GPS, it appears that everything has been slightly shifted. Using one point as an example, in the GPS it is X = 327398 Y = 1836319, but in QGIS it is X = 327401 Y = 1836410.

Finding this strange, I tried using the 'Project' tool in ArcMap to convert the GPX from WGS84 to NAD27 UTM16N. Using this method, my points have the same coordinates as in the GPS.

Does anyone know how to solve this issue in QGIS, or have any insight as to why the coordinates for points are not matching up?


texts: Optimizing PostGIS ST_Buffer to run in near real time

texts:
I have two tables in PostgreSQL 10.1 (with PostGIS 2.3).
The first table contains 7.5 million points (WGS84 coordinates). The other one contains linestrings (protected areas near the roads, 27300 linestrings, wgs84).
I would like to write a query which returns if the selected point (for example where gid=123) is near an area or not (boolean type). It is a very important criteria that the query should run very quickly (1 second or better).

I wrote this query, it use ST_Buffer with 2 meters. It works good, but slowly (3 sec):

SELECT p.gid,
       CASE
          WHEN ST_Intersects(ST_Buffer((p.geom::geography),2),pr.geom::geography)='true' THEN TRUE
          ELSE FALSE
       END
FROM points p
JOIN protected_lines pr
ON ST_Intersects(ST_Buffer((p.geom::geography),2),pr.geom::geography)
WHERE p.gid=3905591;

I set GIST index (and cluster) to both tables.

How can I speed up my query to run in near real-time?

I tried this with polygons (about 1500 pieces) without buffer and ran quickly (better than 100 ms).


texts: Quickly changing outline color of polygon in ArcGIS Desktop without opening Symbol Selector?

texts:
In ArcGIS Desktop 10, I know it's possible to quick change the fill color of a polygon by right-clicking the symbol box of any layer and just pick a simple fill color.

Is there a way to change the outline without open the "Symbol Selector" window?

Sei que é possível alterar a cor de um preenchimento clicando com o botão direito na caixa do símbolo de qualquer layer e aí basta escolher uma com simples.

Tem como alterar o contorno do polígono sem abrir a janela "Symbol Selector"?


texts: Output shapefile has no records after joining shapefile with table and exporting result?

texts:
I have successfully joined a .xls Excel spreadsheet to a shapefile. The join is successful (the number of remaining records seems accurate), although it takes a really long time (like 20 minutes). For reference it's a shapefile of 140,000 parcels and a .xls containing about 18,000 unique records.

When I go to export this newly joined layer so that it is a permanent layer (Data -> Export Data), it works almost immediately and asks me if I want to add the new layer to my map. When I do however, there are 0 records in this new shapefile, but the shapefile has been created and shows up in my computer and in the TOC. There is just no information in it.

The .xls has no spaces/special characters in the column names. Can you think of anything else that may cause this undesired result?


texts: GeoPandas ParseException: Invalid HEX char error

texts:
I am trying to read a geospatial table from a PostGIS database on PostgreSQL.
Using the funcion

sql = '''SELECT ST_AsText(geom) AS geom  FROM my_table;'''
df = gpd.GeoDataFrame.from_postgis(sql, con, geom_col='geom' )

I got back this error

ParseException: Invalid HEX char
Traceback (most recent call last):
  File "nuts3_grid_intersection.py", line 35, in <module>
    df = gpd.GeoDataFrame.from_postgis(sql, con, geom_col='geom' )
  File "/usr/local/lib/python3.6/dist-packages/geopandas/geodataframe.py", line 274, in from_postgis
    index_col=index_col, coerce_float=coerce_float, params=params)
  File "/usr/local/lib/python3.6/dist-packages/geopandas/io/sql.py", line 55, in read_postgis
    geoms = df[geom_col].apply(load_geom)
  File "/usr/local/lib/python3.6/dist-packages/pandas/core/series.py", line 3194, in apply
    mapped = lib.map_infer(values, f, convert=convert_dtype)
  File "pandas/_libs/src/inference.pyx", line 1472, in pandas._libs.lib.map_infer
  File "/usr/local/lib/python3.6/dist-packages/geopandas/io/sql.py", line 54, in load_geom
    return shapely.wkb.loads(str(x), hex=hex_encoded)
  File "/usr/local/lib/python3.6/dist-packages/shapely/wkb.py", line 14, in loads
    return reader.read_hex(data)
  File "/usr/local/lib/python3.6/dist-packages/shapely/geos.py", line 409, in read_hex
    "Could not create geometry because of errors "
shapely.errors.WKBReadingError: Could not create geometry because of errors while reading input.

I am not able to find a solution for it.
I have geom data on a WKB format. I had the same issue also if I convert them using the functions ST_AsBinary(geom) and  ST_GeomFromWKB("geom",SRID).

I run the PostgiGIS code:

SELECT ST_GeomFromWKB(decode('geom ', 'hex'), 4258) FROM my_table;

I got that output

ERROR:  invalid hexadecimal digit: "g"
********** Error **********

ERROR: invalid hexadecimal digit: "g"
SQL state: 22023

How can postgis read the geometry data even if the HEx standard on the WKB is not correct?

How could I solve this?


texts: Is there a difference in using Haversine distance vs re-projecting coordinate frame?

texts:
I have a GeoTIFF raster file in the EPSG 4326 (WGS 84) geographic coordinate system.

I want to crop a 20 km x 20 km square around a particular (lat, long) coordinate, and save that as a new raster.

But since the units for my image are in lat long degrees, but I want a box in terms of kilometers, is the best way to do this to use the Haversine formula to convert my (lat, long) degrees to coordinates in meters, compute the bounding box in meters, and then invert the Haversine computation to get back the box in degrees?

OR do I just project the entire image into a projected coordinate system, so that I can then easily find the bounding box in meters, and then re-project all the bounding box coords back into degrees in my original geographic coordinate system?

Is there any difference in accuracy when doing this? I understand that using Haversine may not be that accurate, but it seems to me that projecting would also lose some accuracy?

(I'm doing this in Python with rasterio or GDAL if that makes a difference.)

My data is around lon: -120.011597902, lat: 39.54719334000001
Resolution (cell size): 0.0035 deg/pixel
- I think this is around 300 meters/pixel (not sure).


texts:
I have a dataset which I believe is generated for ArcGIS .It includes a .tif, .tfw, .tif.aux, .tif.ovr and a .tif.vat.dbf files.

I want to view this in R like other users view it with ArcGIS.

Is there a way?

I searched online, but I couldn't find a systematic answer. Some one says using raster but didn't give any code, and some one says there is a tool package, but no any further detail.


texts: Checking if network has connection errors

texts:
I have given a network of a city in dwg.
I made a dxf that only contains the network itself, and created a layer for every type of pipe (like 25KPE, 80KPE etc...).

How do I check if the network is drawn correctly?

I try to find every error I can, but im afraid the services pipes are either just a little bit shorter, or longer, and not connecting to the main pipe properly.

I'll have to use this network in EPANET later.


texts: Clip layer with selected object from another layer in QGIS

texts:
In ArcMap, I can clip a layer with selected objects from another layer by using the editor toolbar.

I want to do the same operation in QGIS. I want to create holes in the landuse layer by selecting objects from the lake layer. Is it possible?

I can't find a similar option in QGIS 2.8.1.


texts: WMS GetFeatureInfo implementation over Google Maps

texts:
I am using the Google SDK for my base maps. I can see my layers from GeoServer overlayed on Google Maps but I am stuck at the GetFeatureInfo request.

(void) mapView: (GMSMapView * ) mapView1 didTapAtCoordinate: (CLLocationCoordinate2D) coordinate {
    CLLocationCoordinate2D position = CLLocationCoordinate2DMake(coordinate.latitude, coordinate.longitude);
    NSLog(@"Latitude= %f lonitude = %f", coordinate.latitude, coordinate.longitude);
    GMSMarker * marker = [GMSMarker markerWithPosition: position];
    point = [mapView1.projection pointForCoordinate: marker.position];
    //ankita.dutta@hardcastlegis.com
    GMSCoordinateBounds * bounds = [
        [GMSCoordinateBounds alloc] initWithCoordinate: coordinate coordinate: coordinate
    ];
    NSLog(@"Bounds.minx = %f,%f,%f,%f", bounds.southWest.longitude, bounds.northEast.latitude, bounds.northEast.longitude, bounds.southWest.latitude);
    NSLog(@"Pointer: %d, %d", [
        [NSNumber numberWithFloat: point.x] intValue
    ], [
        [NSNumber numberWithFloat: point.y] intValue
    ]);
    GMSVisibleRegion region;
    region = mapView.projection.visibleRegion;
    GMSCoordinateBounds * bounds1 = [
        [GMSCoordinateBounds alloc] initWithRegion: region
    ];
    NSLog(@"Bounds.minx = %.15f,%.15f,%.15f,%.15f", bounds1.southWest.longitude, bounds1.northEast.latitude, bounds1.northEast.longitude, bounds1.southWest.latitude);
    [self getInfo: [
        [NSNumber numberWithFloat: point.x] intValue
    ]: [
        [NSNumber numberWithFloat: point.y] intValue
    ]: bounds1.southWest.longitude: bounds1.northEast.latitude: bounds1.northEast.longitude: bounds1.southWest.latitude];
}
(void) getInfo: (int) xValue: (int) yValue: (double) bboxLMin: (double) bboxLMax: (double) bboxRMin: (double) bboxTMax {
    // [self showLoadingView:YES];
    dispatch_async(dispatch_get_global_queue(DISPATCH_QUEUE_PRIORITY_DEFAULT, 0), ^ {
        //2048x1536
        NSMutableURLRequest * request = [
            [NSMutableURLRequest alloc] initWithURL: [NSURL URLWithString: [NSString stringWithFormat: @"http://ipaddress/geoserver/wms?service=WMS&version=1.1.1&srs=EPSG:4326&bbox=%.15f,%.15f,%.15f,%.15f&styles=&&buffer=40&info_format=application/json&request=GetFeatureInfo&layers=gidc:plots_category_wise&query_layers=gidc:plots_category_wise&width=%d&height=%d&x=%d&y=%d", bboxLMin, bboxLMax, bboxRMin, bboxTMax, 1536, 2048, xValue , yValue ]]
        ];
        NSURLSession * session = [NSURLSession sharedSession];
        [
            [session dataTaskWithRequest: request completionHandler: ^ (NSData * data, NSURLResponse * response, NSError * error) {
                if (data != NULL) {
                    dispatch_async(dispatch_get_main_queue(), ^ {
                        //[self showLoadingView:NO];
                        NSMutableDictionary * listDict = [NSJSONSerialization JSONObjectWithData: data options: 0 error: nil];
                        NSMutableArray * dataArr = [
                            [listDict valueForKey: @"features"] valueForKey: @"properties"
                        ];
                        [self response: dataArr];
                    });
                }
            }] resume
        ];
    });
}


texts: GDAL doesn't recognize gdal.ViewshedGenerate (Python API)

texts:
GDAL has a function called gdal_viewshed. I checked on the GDAL's Python API to see if there is a function I can call in a Python script and I found gdal.ViewshedGenerate.
https://gdal.org/python/index.html
However, when I run a script calling this function I get the following message:
AttributeError: module 'osgeo.gdal' has no attribute 'ViewshedGenerate'

I tried with the latest GDAL version (3.0.2) and version 2.3.3 and I always get this message.
The problem is not the code but the fact that GDAL doesn't recognize this function. Does anybody know why I can't use this function although it is on the Python API documentation?
I couldn't even find any information about this function anywhere apart from the mentioned documentation.


texts: QtDesigner shipped with QGIS 3 doesn't show svg files preview in resources

texts:
In the QtDesigner shipped with QGIS 3 svg's and jpg's are not shown in the Resources window nor as icons on widgets:

QtDesigner from QGIS 3.2:

svg's example

jpg's example

QtDesigner from QGIS 2.18

svg's example

jpg's example

In a QGIS 3 plugin that uses the .ui file the icons on any widget display correctly, so the result is ok, but working with the QtDesigner that comes with QGIS 3 is hardly possible in my case.

Is there a way to make it show my svg's and jpg's as usual, and if not, is it critical or unsafe to use QGIS 2's QtDesigner for QGIS 3 plugin guis?


texts: How to find Dryest Pixel of an ImageCollection in Google Earth Engine? (Dryest pixel composite)

texts:
I want to create a composite with the dryes pixels of an ImageCollection. I know how to make a greenest pixel composite, but I have no clue how to make a dryest pixel composite.

Is there a function similar to ee.ImageCollection.qualityMosaic(qualityBand), but for the dryest pixel?

Code for greenest:

var collection = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA");
var filtered = collection.filterDate('2015-04-01', '2015-07-01');

// Greenest Pixel Composite
function addNDVI(input) {
 var ndvi = input.normalizedDifference(['B5', 'B4']);
 return input.addBands(ndvi);
}

var with_ndvi = filtered.map(addNDVI);
var greenest = with_ndvi.qualityMosaic('nd');
Map.addLayer(greenest, {max:0.3, bands:["B4", "B3", "B2"]}, 'RGB');


texts: Selecting multiple partial field contents using ArcGIS Select By Attribute?

texts:
I can't seem to get the hang of ArcGIS Select By Attributes SQL. I'm trying to "select all that contain Street and Road and etc." Working with polyline SHP.

This works:

"NAME" LIKE ('%Street%')

This doesn't:

"NAME" LIKE ('%Street%', '%Road%')


texts: Custom Certificate Configuration error message

texts:
When I run QGIS v. 2.18.15 on my PC I get this error message when the program first starts up:

Unable To Get Local Issuer Certificate: The issuer certificate of a locally looked up certificate could not be found

Is there any way I can fix this issue, because currently I am not able to install any new plugins. I've tried googling the error message but I am not able to find any helpful solutions!


texts:
First of all, it should be understood that I am a complete novice to javascript and interactive maps.

I'm working on building a web map that will show a time lapse of all of the worlds nuclear detonations from the years 1945 to 1998. Simply put, there is a map and a play button with a slider next to it. When you hit play, markers will begin to be added to the map as the timeline progresses (it takes about a minute in total for the timeline to go from 1945 to 1998). Setting the timeline slider to a specific date will have the map display a marker for every nuclear detonation that had happened up to that date.

At the moment, the markers are being added to the map properly after hitting play, but my problem comes when I drag the timeline slider backwards. The map will always keep the markers of the previously set date up on the map. Clicking the play button again will continue the timeline and make duplicates of the previously added markers.

I am thinking that I might have to do some sort of if-then statement that says something like if the date = x, remove all markers past date x. Any pointers?

I'm using quite a few plugins / js-libraries. Here are the ones that show up in my script:

Leaflet - for mapping

Leaflet.timeline - timeline plugin (https://github.com/skeate/Leaflet.timeline)

Leaflet.markercluster - for clustering together blobs of markers

Moment - for formatting dates found in my dataset to be used in the timeline

jquery - for jquery stuff

And here's the actual timeline portion of my script:

function sipriReportTimeline(data){
var getInterval = function(feature) {
    return {
        start: feature.properties.date_long_moment,
        end: moment("1998-05-30"),
    };
};
var timelineControl = L.timelineSliderControl({
    formatOutput: function(date){
        return new Date(date).toLocaleDateString();
    },
    duration: 100000,
});
var timeline = L.timeline(data, {
    getInterval: getInterval,
    pointToLayer: sipriReportPointToLayer,
});
timelineControl.addTo(myMap);
timelineControl.addTimelines(timeline);
timeline.addTo(myMap);
};


texts: Qgis - unwanted offset when using raster interpolation plugin

texts:
I'm using QGis and I have a set of points with a depth attribute. I'm trying to make an interpolation using the Raster interpolation plugin (Raster -> interpolation -> interpolation). But when I do so I get a misplacement on the produced raster. It's positioned too far north. Do you have any ideas of what may be causing this?

To make it clear what I am doing and what my problem is I've created a short video where I'm showing the issue:
http://jimutt.com/interpolationProblem.swf

As you can see in the video I'm using the same projection in both the viewport and for the point layer.

Please tell me if you can't display the video.


texts:
<Rule>
    <PointSymbolizer>
     <Graphic>
       <ExternalGraphic>
         <OnlineResource
           xlink:type="simple"
           xlink:href="glass.png" />
         <Format>image/png</Format>
       </ExternalGraphic>
       <Size>28</Size>
     </Graphic>
   </PointSymbolizer>
</Rule>

I want to use external graphics on my layers. I styled with geoserver admin panel. But external graphics are not working. I put my image into /var/www/geoserver/data/styles folder.

I tried some href options but not worked.

xlink:href="glass.png"
xlink:href="file:///glass.png"
xlink:href="/var/www/geoserver/glass.png"

I think another issue.


texts:
Three and a half years ago a user asked if it is possible to add bookmarks to the print composer in QGIS. The link is here.
The answer was that someone posed a feature request for this feature.

I was wondering if this feature is now available? Perhabs it came with QGIS 3? I can't find it though but maybe I just don't see it.


texts: Why a laz file warns: 'return number' greater than 'number of returns' in R?

texts:
Reading a laz file with the lidR package in R I have the following warning. What does it means?

las = readLAS("chmp15.laz", filter = "-drop_z_below 0")
#> Warning message:
#> Invalid data: 8895 points with a 'return number' greater than the 'number of returns'.


texts: Problem deleting parts / frame of polygons / shape

texts: I generated a polygon-layer from a raster dataset. the extensions of the raster are now shown as a frame. When I delete this and save it is there again (pls see attached images): 1st image how it looks before deleting frame and after saving, 2nd image how it should look like. Can anyone help me on this???


texts: Can I color different layers with the same pseudocolor color map?

texts:
I have multiple DEM layers in my map (downloaded from atlas.louisiana.gov/lidar), but I cannot seem to get the colors to line up.  You can see in the attached image that the colors are not continuous between the layers, but they should be.  (I've checked the values, and they are).  I've tried several ways around this, but none of them seem to work.  This data is rendered as Single Band Gray.

1)  right clicking on one layer, selecting copy style, and then pasting it into the other layer.

2)  saving the style from one layer and loading it into the other.

3)  changing the color map to grayscle, inputting min and max, and then reverting back to pseudocolor.  (that one seemed silly)

It seems like the problem is that I can only do a pseudocolor with standard deviations, and not min/max values.  Is that correct?  There seems like there should be some way around that...

Any ideas?


texts:
I'm trying to make multi-layer vector tile using SQL query. I could able to make Mvt tile, but that was only for single layer. I got binary file from that query.

How can I make multi-layer Mvt vector tile? I need the same binary file.


texts: Preventing label text from overlapping in ArcGIS for Desktop?

texts:
Using ArcGIS 10.1 Desktop of Home Use (ArcInfo license I believe), I'm getting results like this using the included Maplex labelling engine:

Placement is curved. [Try horizontal position first box is unchecked.]

Strategy order for fitting is:

Stack Label - default settings
Reduce label in size - Font is 12 pt Arial, lower limit set to 1 pts, step interval is 1.1 pts.  Font Width compression is at default (90% lower limit & 5% step interval).
Abbreviate Label - not enabled
Overrun feature - not enabled
Compress label in width - not sure if this enabled or not.

The region "Steamboat - Hilltop 16" is properly labeled.  As long as the entirety of that text is within the polygon in question I don't care how it might be stacked.

The region "Steamboat - Mountain 17" (directly below Hilltop 16 and to the left) is incorrectly labeled.  It has written the "- 17" underneath the Mountain and is unreadable.

This didn't happen using Maplex in ArcGIS 10.  Is there anyway to stop this behavior?


texts: Can't find protected_areas.shp in qgis_sample_data or GIS database

texts:
I have just installed QGIS 3.10 and have unzipped the two sample files named in heading cannot find the files named in:

https://docs.qgis.org/3.4/en/docs/training_manual/introduction/preparation.html

and have searched online for protected_areas.shp and exercise_data/shapefile/protected_areas.shp ?


texts: Display attribute values over time in a graph as a widget or pop up in ArcGIS Online

texts: I'm trying to find a way to display weekly settlement monitoring data from a set of points as a line graph in ArcGIS Online.  The charts functions that I've played around with require a statistical operator which is not helpful to the data I am trying to display.  Ideally, each week the elevations of a series of points would be joined to the existing hosted shapefile and the graph would update automatically. Column headers would represent dates and the attribute values would be the elevations.  Am I missing an available function/tool or is this not a possibility in ArcGIS Online yet? Am I going about this in the wrong way?


texts: XYZ Clustering in Python and numpy arrays

texts:
I am looking for a method of clustering XYZ data in Python that works on numpy arrays.  I know of a few sources, such as clusterpy and Pysal but have had little success with them as they seem to still be in the development stage and do not support numpy arrays.

EDIT:

My numpy arrays are of XY points with a Z column of a value at the point such as soil moisture or plant height.  I want to spatially cluster the Z values into similar regions.


texts: Open Source Address Correction / Parser with Fuzzy Matching

texts:
Here is a a little bit of a detailed question related to address parsing/geocoding which I feel should be interesting to many users.

So, essentially I am curious to know if anyone has had any experience installing, building or extending a opensource geocoding and/or address correction tool.

I am aware of geocoder:US 2.0 initiatives which I think are maintained by geocommons but I am unsure if there are better alternatives, other opensource tools, if their system can be effectively extended or if there are any developments I might not be aware of.

My goals are as follows:

I need a highly accurate tool which is capable of automatically parsing out and/or standardizing location data inputted by users from a single input field all in real-time and with the highest volume possible.
Input data would be one or more of the following address components: zipcode, county, city, street, address, state.
Input data also needs to be able to lookup from our custom geonames database. For example he may enter the name of a neighborhood or non USPS location name which naturally are not standard address variables.

Given these goals I am well aware of the fact that when given a single form field to conduct such a lookup each user will enter his data in different formats while the other factor generally fall into misspellings.

Besides utilizing the census database as the core for the valid addresses/ranges (all which I believe Geocoder:US does, I believe some type of ability to define known "aliases" would be ideal for known misspellings of street names. The same goes for things such as a user entering Ave compared to Ave. compared to Avenue. Don't think such alias capabilities are fully possible with the Geocoder:US tool.

While the above elements may indeed solve the majority of issues I think some type of effective fuzzy matching needs to exist when the input can't be matches to high enough %age.

If input data can effectively be parsed out into individual elements based off some assumed rules and then utilizing a type of "match score" component to fuzzy match any unmatched elements would have to be based on those elements which were already "matched" with a high degree.

For example: I am going to assume for geocoding to be as effective as possible we need to extract individual data elements from the input field first in an attempt to narrow down the "area" the user is trying to find results for. In my view this means that a 5 digit number could be assumed to be a zipcode, if there is another element such as a city name that matches the zipcode the assumption that we have the "area" correct... Next we use the remaining data to try to find a full, partial or fuzzy match, score and list possible results.

In any case - I would greatly appreciate if anyone could provide some advise here along with any advise, performance stats or upcoming developments they are aware of which might adjust my direction (such as the use of postgis 2.0 as a means for enhanced matching capabilities)


texts: How to merge two Polygons in OpenLayers?

texts:
I have two polygons.

polygon1 = new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Polygon([linearRing1]))
polygon2 = new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Polygon([linearRing2]))

How can I merge the two polygons into one using OpenLayers?


texts: Multiple tabbed panels in Google Earth Engine app?

texts:
How can I add multiple tabbed panels to a Google Earth Engine app?
I don't see an option for this in the docs.
My app has two components, data visualisation, and some options to download points for ground truthing. All the components together means users would have to scroll up and down quite a bit in the panel, but having multiple tabbed panels would solve this problem.


texts: How to use wglUseFontBitmaps in ArcGlobe?

texts:
Does anyone have any experience using, wglUseFontBitmaps in ArcGlobe.

It seems no matter what I do, the text doesn't show up on the globe.

Here's my code.  This is the same way that I've plotted points, drawn polygons, polylines, etc on the globe.  Any help would be great.

double outX = 0;
double outY = 0;
double outZ = 0;

myGlobeViewUtil.GeographictoGeocentric(0 ,0 ,0 , out outX , out outY , out outZ);

Font font = new Font("Courier New" , 32.0f , FontStyle.Bold);
uint dbase = GL.genLists(256);
wglUseFontBitmaps(this.Handle, 32, 256 , dbase);

GL.glLoadIdentity();
GL.glTranslatef(0.0f , 0.0f , -1.0f);

GL.glColor3f(1.0f , 0.0f , 0.0f);
GL.glRasterPos3d(geoX , geoY , geoZ);
GLPrint("TEST");

private void GLPrint(string inText)
{
  GL.glPushAttrib(GL.GL_LIST_BIT);
  GL.glListBase(dbase -32);
  GL.glCallLists(text.length , GL.GL_UNSIGNED_SHORT , text);
  GL.glPopattrib();
}


texts: Set one output band's no data value with gdal_merge.py

texts:
I'm trying to stack two images with gdal_merge.py -separate so that I can use the NoData values in the second band to mask the first band (was going to use gdal_translate -b 1 -mask 2). The first image (or band 1 after the merge) has no NoData values, while in the second image NoData = 0. However when I use gdal_merge to stack the two images, the NoData for the second image/band is removed and set back to 0.

There is a "-a_nodata output_nodata_value" option in gdal_merge, but I don't want to set the first band's NoData values, just the second. Is this possible?

Or is there another way I can use the second image's NoData values to mask out values in the first without stacking the images via merge? That is my ultimate goal. Maybe I can achieve this via gdal_calc?


texts: Converting CSV file of WKT data to shapefile using ogr2ogr?

texts:
This question is related to Converting shapefiles to text (ASCII) files?.

I have a CSV file, with one column, where all the rows correspond to WKT POLYGON()'s:

WKT
POLYGON(...)
POLYGON(...)
...

I'm familiar with how to convert from shapefile to wkt, but I need to go the other way around. How to accomplish this task?

I tried playing around with ogr2ogr's settings/flags but didn't really get anything useful.

I also know I can use QGIS to do this, but it freezes/crashes since the dataset is fairly large.


texts: Help with joining a large DBF file to shapefile

texts:
My overall objective is to create a time line animation depicting when a person is working in a particular region.  The .dbf table consists of over 6000 record entries of shift data, including date start and date stop times along with an assignment ID.  I have approximately 100 different shapes/regions also with assignment ID's.  I would like for that region to turn on or off based of the table shift data and show when someone is working that region.

When I join the table to the shapefile, only the first shift entry gets attached and displayed.  Do I have to have 6000 shapes for each of the 6000 record entries?

There has to be a better way, please help.


texts:
I'm having a bit of a memory leak, or more like a memory sieve really.

I'm working in C++ and I'll illustrate what's happening with one iteration of part of my program.

I have an IGeometryPtr, and for this case it is pointing to an IPolygon with ~700,000 points in it (the issue is present, but less pronounced with smaller geometries). All is good and well until I attempt to extract point or ring information from it. At that point my memory use spikes way up (which is to be expected at this point). The issue is that a chunk of that memory is never released. I'll usually lose about 100 MB on this gigantic geometry. Here's the code:

void ShortLengthSegmentHelperNative::loopPolygonParts(IPolygon4Ptr polygon, void* arg1, void* return1, void (ShortLengthSegmentHelperNative::*f)(IGeometryPtr, void*, void*, bool&))
{
    long extRingCount;
    polygon->get_ExteriorRingCount(&extRingCount);

    //If I return here there is no leak!!!

    IRing** extRings = new IRing*[extRingCount];
    polygon->QueryExteriorRingsEx(extRingCount, extRings);

    //Leak appears to be the same size regardless of releasing the Rings
    for (int i = 0; i < extRingCount; i++)
    {
        extRings[i]->Release();
    }
    delete[] extRings;
    //lose 100 MB when returning here
    return;
    //...stuff to do when this is fixed
}

I should note that this also happens if I retrieve the geometry bag instead of using QueryExteriorRingsEx. The polygon shouldn't need to be released because it's a smart pointer and will release itself.

Does anyone know why this code leaks memory? This is ArcObjects 10.2


texts: I have point layer with text field. There is many empty field. I want to select only filled text. Something like this, with QGIS:


texts:
I have two different files. The first file is from the United State Census. It is a file consisting of the census tracts for 1990. Each polygon has information that labels its census tract, namely a variable named "TRACTBNA".

The second file is from the King County GIS office. It is a shapefile of all the property parcels within King County. Each polygon has a unique PIN (parcel identification number).

Both the parcels and census tracts are polygons. I need each PIN or polygon from the second file to be given a TRACTBNA. I used join by using the join and relates, and choosing "join data from another layer based on spatial location". I also chose the second option.

I exported the data as dbf and I opened it with Excel. I find that not every parcel was given a TRACTBNA, in other words it is empty.

My guess is that some parcels may not lie completely inside a single polygon, and those aren't spatially joined. I don't care if it randomly selects one of the polygons it intersects, as long as all of them receive something.

I am using ArcGIS.


texts: QGIS: RiverGIS: Couldn't connect to river database

texts:
Trying to connect to the PostGIS DB for RiverGIS fails with
Couldn't connect to river database: var@localhost
Please, check you database connection settings!
This happens despite the DB being connected to/loaded (see screenshot below) and already having loaded a layer from the DB into the currently open project which I can edit without problem.
What am I missing respective what steps can I take to debug this?
Edit: Having tried answer 1, adjusting the the user & db-name to be same as in the documentation does not solve the problem.
System: Win10 / QGIS OsGEO4W / 3.18.3
I'm aware of QGIS - error connection of database in RiverGIS but there the DB doesn't seem to have been connected to...
Edit: Adjusted screenshots, having tried answer 1:


texts: QGIS point displacement does not work as intended

texts:
Maybe I made a mistake, but I can't get it working. I have a CSV file with species-sightings and lot of them at the same locations. So to show all species I need to use point displacement. If I choose to show them categorized the outcome is as expected(pic 1), they are all over the map. But if I choose point displacement, all of them are around just one point (pic 2), that seems to be in the middle of the location (it's not the whole map in the pics).

It seems like I need to specify somehow the coordinates again? Or is it maybe a bug?

QGIS-Version: 2.14.3


texts: Field editor in QGIS 3.0 no Text edition button

texts:
I used QGIS 2.18 and now I use QGIS 3.0 under Windows 10 in Spanish versions.

In the field editor in version 2.18 each field had his button "Text edition" but in version 3.0 that button doesn't appear. In the next images you can see what I mean, the first is from version 2.18 and the second is from version 3.0.

Can someone help me? Have they moved that option to another window?

https://imgur.com/JgCOIAH

https://imgur.com/SjTwPbd


texts: Parsing .dxf through a request to the ArcGIS server REST API

texts:
I have data in the DXF file.
Can I parse the data and get response in an json by querying to the ArcGIS server (REST API)?


texts: Which operations must be performed using the Geoprocessing API in ArcGIS Pro?

texts:
I want to migrate .NET application from arcobjects to Arc GIS Pro SDK.
As I understand some operations such as creating tables and feature classes, creating and modifying fields need to be performed using the Geoprocessing API in Arc GIS Pro SDK.
Where can I find full list of such operations or more information about it?


texts: how to output convex_hull polygon of points shapefile

texts:
i have a shapefile of more than 8 points and i have written a code get the convex_hull geometry of of the points and output it as a shapefile using fiona. however the code returns an error saying:

ERROR:TypeError: expected string or bytes-like object

from shapely.geometry import MultiPoint, point, shape
from fiona.crs import from_epsg
outshape = fiona.open(".cache\ca2_points\ca2_points.shp")
with fiona.open('.cache/mew.shp', 'w', 'ESRI Shapefile', outshape.schema.copy(), outshape.crs)   as output:
    for feat in output:
        feat = MultiPoint([shape(point['geometry']) for point in fiona.open(outshape)])
        output.write(feat)


texts:
I'm trying to convert osm.pbf to .shp for some analysis in QGIS. I need to verify id the converted .shp is the same as .pbf.
But the attributes table shows only 99 records by default. How do i find the total number of entries?


texts: Making animated choropleth with ArcGIS Pro

texts:
I want to make an animation that shows the changes in state shape and state population every ten years. My state data is all in one large feature class. I have a polygon for each state every ten years, along with the state's population.

I want the animation to start with something like this:

and end with something like this: I want the animation to cycle through the different years, with maybe a 2 second pause for each year.

I know how to use graduated colors in the symbology to make the choropleth maps, I just don't know how to use the animation tools in ArcGIS Pro to make the animation.

I tried to create key frames of the map using different Definition Queries (ie CensusYear=1900), but the key frames wouldn't change. Do I need to make separate maps for each key frame? What is the best way to do this with ArcGIS Pro?


texts: QGIS 2.14.17 repositories "unavailable" cant install plugins

texts:
I cant access the plugin repositories. If i go to plugins -> Manage and install plugins and settings i get unavailable for QGIS official repositorie with url "http://plugins.qgis.org/plugins/plugins.xml"
How can I fix this?


texts: Geoserver linestring fill based on attribute value

texts:
I'm experimenting with Geoserver. I'm styling a linestring layer with CSS in this way, which works fine:

[userid = 1]  {

stroke: #8C2300;
stroke-width: 3px;
stroke-linecap: round;
z-index: 0, 1;

}

In my data store I have an attribute 'color' with hexadecimal values. Can I use this value for stroke color, in this way:

[userid = 1]  {

stroke: [color];
stroke-width: 3px;
stroke-linecap: round;
z-index: 0, 1;

}


texts:
I'm new to QGIS and I'm using it as a small part of a project I'm developing, so I need some help if possible. I have this layer:

What I need is to define a variable to serve as a perimeter and give some lat/long coordinates and define a circle around that point with perimeter = variable. I don't need this to all the points on the layer, just a few that I will then select. Also, the perimeter should be a constant. The main goal is to count the number of red/yellow/blue "points" within the range of each circle.


texts: Adding xy data to ArcGIS Desktop and converting it to grid that is clippable?

texts:
I have an Excel file with x, y coordinates (lat, lon) and also a value column for my index of interest.
I added the data using

  file->Add data->Add xy data

Then I want to produce a 'grid' of this point data I have added (to produce a nice map of course).
I then used the 'feature to raster' option within the conversion toolbox, and I set the grid size to 0.5 (what I'm interested in). This produces a grid of my data with the values categorised by colour.
Then, I try to clip it using a world map coastline using

  Raster processing->Clip

in the 'Data management'toolbox.
It never seems to work however. All it seems to do is produce a new raster (grayscale) of the data, with nothing clipped out around the coastline boundaries.


texts: point returned from ST_LineLocatePoint returns false in ST_Contains

texts:
I am using postgis's ST_LineLocatePoint to find out the closest point on a LineString to the given Point, and using ST_LineInterpolatePoint to extract a Point from the returned float number.(as referred here)

This is my ST_LineLocatePoint Query:

SELECT ST_AsText(ST_LineInterpolatePoint(foo.the_line,
         ST_LineLocatePoint(foo.the_line,ST_GeomFromText('POINT(12.962315
         77.584841)')))) AS g FROM (SELECT ST_GeomFromText('LINESTRING(12.96145
         77.58408,12.96219 77.58447,12.96302 77.58489,12.96316 77.58496,12.96348
         77.58511)') AS the_line) AS foo;

Output:

                  g
------------------------------------------
 POINT(12.9624389808159 77.5845959902924)

Which exactly lies on the linestring I have passed. Demonstration is displayed here.

But when I check whether this point lies in the same linestring using ST_Contains it always return false, even though the point lies within.

My ST_Contains Query:

    SELECT ST_Contains(ST_GeomFromText('LINESTRING(12.96145 77.58408,12.96219
    77.58447,12.96302 77.58489, 12.96316 77.58496, 12.96348
    77.58511)'),ST_GeomFromText('POINT(12.9624389808159 77.5845959902924)'));

Output:

 st_contains
-------------
 f

I am not getting where I am doing wrong, why I am getting false. Can anyone help me in this.

Postgresql : 9.4
postgis : 2.1


texts: Extracting Factors Causally Related to Landslides from DEMs (SRTM)

texts:
I have the SRTM dataset, and I've managed to write some python code to show elevation and slope heatmaps. The code runs fine, and the elevation data is reliable (I checked some points from the web), but the slopes that I'm getting seem to be really off. I checked my code and it doesn't seem like there's any logical error there.
I'm using python's richdem library for Terrain Calculations: https://github.com/r-barnes/richdem
And my code for calculating slope is here: http://pastebin.com/1PxWMJ4A.
Here's an example elevation map on the left and the corresponding slope map on the right (degrees). (they're a bit small so feel free to open in a new tab).

The region I've chosen is super flat (with the exception of the river and the hills up north).
Why is the resulting slope map so noisy?
Is it supposed to be like this?
If I need to take a measure of the slope around this area (say, to assess landslide hazard), how am I supposed to do this without having my metric be subject to all the wonkiness in the slope map?
Additionally, here's a global slope map that I generated (again, in degrees). So much of the earth's surface has a slope of >50 degrees!
What is going wrong?

For example, the following image is a place near the Himalayas where landslides are common.

How can I get any useful information if I'm being blinded by slopes >= 80 degrees all the time?
Am I just calculating the slope incorrectly?
I did calculate the rise-run slope for the same area. (see below) It looks like it's better, and I'm willing to use it, but I want to understand why the slope in degrees is being wonky.


texts: Using Make Feature Layer tool in ModelBuilder?

texts:
I am having a slight problem with using the Make Feature Layer tool in ModelBuilder.

Basically I have a file geodatabase of points. I am trying to make a model to pull out points according to a reference field value. I could just make the expression for the tool a variable but unfortunately the users who are going to use it are to dense to create a simple SQL statement.

I have created a model String variable (VRefNo) to take the reference number (Its alphanumeric hence the string. The column in the geodatabase is also a string) and tried to parse it in the usual way to the expression:

"RefNo" = '%VRefNo%'

and all the other variables are hard coded. I have made the variable a model parameter and checked the add to display check box for the output and made sure nothing is intermediate (There is only one tool in the model and thats it).

The model runs ok with no errors but I am not getting an output.

Any Ideas?

I can get the select tool to work in ModelBuilder when run from the toolbox only when I make the output a model parametre. That's all nice but I only want to save the results temporarily to memory (in_memory) and would prefer the user not having the option to change this.

I would make a python script to do it but all the clients dont have python installed on their machines.


texts:
I have a polygon layer. I drew a boxplot calculated from field area in the layer.
How can I select the features whose area value is within interquartile range (25%-75%)?


texts: Assigning common field value to segments between two points using ArcGIS for Desktop?

texts:
I am using ArcGIS 10.3 for Desktop.

I have a shapefile with polylines connecting points. Instead of having one single polyline connecting two points, there exist a number of consecutive segments connected to each other (see image below).

These segments, however, have no field values in common, so if I wanted to use Dissolve to merge all of them together it would not work the way I expect.

Is there a way to assign a common field value to all the segments between two consecutive points? This way, Dissolve could be used effectively. I am open to any kind of solution, and would prefer a Python-based one, if available.

Visual example: all segments in between point A and point B should get a common field value, so that they can be later merged. This should apply to all segments in between any two points.

This question is not a duplicate of ArcGIS 10.3: How to merge lines divided by junction points?, which I myself have asked. Rather, this question is a precursor to it, if you like.

Screenshot of the real situation. Red dots are the equivalent of Points A and B, while the segments are recognizable by the junction points in green (there is a segment in between each pair of consecutive green dots). For each red dot, I want one single incoming and one single outgoing line.


texts: Displaying ArcScene animation (*.asa) in Powerpoint?

texts:
I just used ArcScene to make an animation clip, whose format is .asa.

Is there a way for me to put it into PowerPoint for the sake of presentation?


texts:
I have developd a map with TOC, Navigation Toolbar and some other query options using Javasctipt API in .aspx page. When I run it individually, it is working perfectly, but I want to embed it into another .net application.

In .NET application I have a page with one dropdown with various options. Upon selection of option from dropdown I need to show the above said map.

Iam getting below error.

**Message: Could not load class 'dijit.layout.BorderContainer
Line: 15
Char: 133114
Code: 0

I have hosted Javascript API in my local IIS.

Can anyone help me how to resolve the above issue.

Thanks,
Uday