texts: Convert several xyz data into a homogenous raster data with GRASS?

texts:
I am new to GRASS and Python and I have to convert a huge amount of *.xyz data into raster data using a batch process in GRASS.
So far it works good, but when I put the created *.tif data into QGIS it looks like some kind of puzzle because of the different colors.
Does anyone if there is a command for GRASS to get the DEM homogenous? I've already tried it with r.rescale.eq but that didn't solve the problem.

My code now looks like this

call r.in.xyz --overwrite -s -g input=%1 output=%2 > tmpRegion
call SET /p myregion= < tmpRegion
call DEL tmpRegion
call g.region %myregion% res=1
call r.in.xyz --overwrite input=%1 output=%2 separator="|"
call r.out.gdal input=%2 type=Float32 output=%2.tif

and this is my result


texts:
How do I invert a shapefile with GDAL? Basically, I have a global water body shapefile. I'd like to ask the question "does this filter region intersect land?" The most obvious way to do this is to invert (or complement) the shapefile and then use the OGRLayer filter functionality.

Unfortunately I can't just go an grab an arbitrary land shapefile -- I have to use this water body shapefile specifically.

EDIT:

I have a water body shapefile:

I would like to invert it (in other words, get it's complement or reverse the polygon winding direction) in to a land mass shapefile:

(Colors aren't important here except for illustrative purposes.)


texts: Are there any spatial interpolation methods that take into account inequality constraints?

texts:
For example if I'm interpolating a water table map from groundwater levels measured in wells and I have a location where a well has gone dry, I don't have a measurement of the depth to water, but I do know that the depth to water is greater than the depth of the well. I would like to have an interpolation algorithm that takes this kind of information into account, so that the value at the dry well location obeys the inequality constraint depth to water > depth of well.

Another example might be if you were interpolating some water-quality or air-quality constituent and had non-detects, you might want to apply the constraint that the value is less than the minimum detection level at those non-detect locations.


texts: QGIS Graphic Modeler not including all GRASS algorithms - How to include v.out.ogr?

texts:
I'm using QGIS 2.18,12 with the GRASS module 7.2.1 on Windows. When building a model in QGIS , I am not able to see all the algorithms possible under thr GRASS module, that are indeed available within QGIS, but simply not in the graphic modeler. Specifically, I need to use v.out.ogr (screenshot attached).

How can I include this in the modeler?

It is vital that I have v.out.ogr built-in the modeler, and not select the function via  Plugins > GRASS, and then from the GRASS Tools: File management> export from GRASS> export vector from GRASS.


texts: SQL-query should return GeoJSON FeatureCollection

texts:
I'm trying to return a GeoJSON FeatureCollection object from a SQL-query with PostgreSQL 9.3 and PostGIS 2.1.

This is my query so far:

SELECT row_to_json(fc) AS p
 FROM (SELECT 'FeatureCollection' AS type, array_to_json(array_agg(f)) AS features
 FROM (SELECT 'Feature' AS type, ST_AsGeoJSON("wkb_geometry")::json AS geometry
 FROM "4547cd00-b570-4504-9796-cf2f61dbebe4") AS f) AS fc;

I'm getting the error:

ERROR:  parse error - invalid geometry

However, if I just run this query:

SELECT ST_AsGeoJSON("wkb_geometry")::json AS geometry
 FROM "4547cd00-b570-4504-9796-cf2f61dbebe4";

It works and print out a nice geometry.
If I ignore the ST_AsGeoJSON it also prints out a result, but not something I can use as a GeoJSON FeatureCollection, as I need the type and coordinates seperated into JSON.


texts: Generating overlapping polygons from lines output from Table To Ellipse tool?

texts:
I have a feature class of line-based ellipses generated by the Table to Ellipse tool that looks like this:
What I want is polygons of the individual ellipses, however tools like Feature to Polygon will create a polygon for every intersection.

My solution to this right now is a Search Cursor in Python that makes a layer from each line ellipse and then runs the Feature to Polygon tool on that layer. This is extremely slow and means I have a feature dataset with hundreds (soon to be thousands) of individual polygon feature classes in it instead of one file.

I saw from this question that there is a Python package, Shapely, that might run much faster than the arcpy tool in a loop. I can't believe, however, that there isn't some way to do this in ArcGIS. Can anyone think of a better solution?


texts:
With QGIS 3.x when I save any project in a qgz format and I try to open it after, they are always with no one layers in the TOC (and canvas...). Every layer is in the same location in the computer...

Anyone has any tip about what could be happening?

Sorry, I wrote the post too quickly...

I am not getting any warning or error at any time
The data is stored locally in the computer hard disk.
The work flow is just:

Add some layer (shp / geopackage / raster)
Save the project
Close QGIS
Open QGIS
Select the project

The OS is Windows 10 (I have tried the same flow in a computer with linux and everything is ok).


texts:
I am very new to GIS in R.

I am trying to split every line in a polyline shapefile into segments of three. My original file has 542 lines, so after the transformation I should have the same shapefile but with 1626 lines. I tried the ArcGis split command but my license only allows me to split one line at the time.


texts: How can I handle NULL values in a color ramp?

texts:
I'm mapping areas in QGIS 2.12 based on average pricing, using area polygons (PostGIS) joined with a CSV data table that doesn't have values for areas that didn't have any sales.

I'd like to color the areas that don't have matching "area" values in my data table, but also use a green to red ramp on the areas that do have pricing.

My data is set up like this:

area, avgprice, sales
A, 50000, 100
B, 100000, 50
C, NULL, NULL
D, 250000,  20
etc.

I was able to set up a zero value with white as the coloring, and 5 other colors encompassing the total average price range, but I can't seem to Classify NULL as a value to use in the color ramp.

Maybe there is a way to use a categorized style, but I don't if there is a way to set value ranges instead of distinct values for that methodology.


texts: In ArcPy, how can I save a clipped raster dataset to memory (not to file)?

texts:
I've been trying to clip a larger raster dataset by a polygon extent, where the resulting clipped raster dataset is saved in-memory. So far, my attempts produce vague errors.

>>> arcpy.Clip_management(in_raster="/path/to/raster_dataset.tif", rectangle="899329.78438381 1065232.74648145 912509.334774774 1080416.77223298", out_raster="in_memory/clipped", in_template_dataset="/path/to/boundary.shp", clipping_geometry="ClippingGeometry")
Runtime error <class 'arcgisscripting.ExecuteError'>: ERROR 999999: Error executing function. Failed to execute (Clip).

Is there a way I can store a raster dataset in-memory?

Documentation for the VB.NET API seem to indicate it is possible beginning in Arc 9.2, and I'm using Arc 10 on my workstation.

  ArcGIS 9.2 allows to create an in-memory RasterDataset where the pixels are stored in memory, this type of raster dataset has the charastics of a file RasterDataset.


texts: Estimating building and vegetation volume with unclassified LiDAR datasets?

texts:
I'm trying to determine the ratio of building to vegetation change (i.e. how has the volume of built structures changed in proportion to the volume of vegetation) in a city over a period of five years. The LAS files for the first time period only have ground, noise, and unassigned class codes; the second period has building, ground, noise, rail, water, and unassigned.

My main difficulty has been creating a building raster for the first time period, which does not have class codes for buildings. Secondly, I need to determine the amount of vegetation in the region for both time periods.

I thought it would be a pretty simple process (perhaps it is), but it has proven rather challenging as I am somewhat new to GIS and Arcmap (10.6). This is for a LiDAR class and the area of interest is only a few sq. km.

Does anyone have any suggestions for how to proceed?


texts:
I've never seen this term, and I'm wondering what it means, with regard to ST_Polygonize

  NOTE: Input linework must be correctly noded for this function to work properly


texts: r raster package rasterize failing:no slots

texts:
I'm having trouble rasterizing a shapefile using the raster package in R.

shp<-shapefile(ZoneShape);
lcRas<-raster(lcRaster);

r<-raster(ncol=ncol(lcRas), nrow=nrow(lcRas), crs=CRS);
res(r)<-res(lcRas);
extent(r)<-extent(lcRas);

>r

class       : RasterLayer
dimensions  : 22610, 27959, 632152990  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 554739, 582698, 3837197, 3859807  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0

shp$GID<-1:nrow(shp);

> shp

class       : SpatialPolygonsDataFrame
nfeatures   : 1
extent      : 554838, 582597.6, 3837297, 3859707  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
nvariables  : 14
names       : SP_ID,              NAME, Shape_Leng, Shape_Area, GID

>zoneRas<-rasterize(shp, r, "GID")

The following error is returned:

trying to get slot "coords" from an object of a basic class ("NULL") with no slots.

Can anyone see what I'm missing up here?


texts: Reduce image collection to get annual monthly sum precipitation

texts:
I'm using GRIDMET image collection in Google Earth Engine, which has daily climate data from ~1980 to present. I would like to create an image collection of monthly (May-September: each month) sum precipitation for each year (1986-2016). I've been trying to create a function to do so, but can't figure it out.

https://code.earthengine.google.com/ceba35ee52955ec5e174e0293cebb7b8


texts:
I am running into an ArcGIS 10.1 problem with two data sets matching up, and although it seems like a coordinate system problem, it has me stumped.

One data set is a GPS track (Garmin Etrex Legend 30) and the coordinate system is set to WGS 1984. It shows what I believe to be the correct decimal degrees in the bottom right corner of the screen. This is a geographic area in Africa (Malawi) that I'm not familiar with so I'm not sure if it's the correct decimal degrees, but when I move the cursor around there isn't too much variation, suggesting it's correct.

The other file was obtained from a state agency and I've used the "Define Projection" tool to use the same coordinate system as the GPS file. It's decimal degrees are very different and seem incorrect due to their discrepancy from the other data set.

When I open the properties for both files and look at the coordinate system they appear the same. Here is the properties for the GPS file

and here is the one for the State's contour file.


texts:
Could someone help me with determining how to export to PDF using QGIS 2.18.13 and my PDF viewer/editor is Foxit?

When exporting to PDFs, it is now blank. I have exported before with no problems. I have tried to do a jpeg, same thing, blank. Exporting to SVG, not blank but text is badly distorted/large in size. I have take text boxes out, still blank pdf. I have updated QGIS, I was using 2.18.1 at first, but I started getting blank PDFs before the update also.


texts: ArcObject10.1：Memory Leak with ILasDataset.Export()?

texts:
I  have created a function for exporting to the las in the geometry range by the ILasDataset.Export method.now every time the function when it is called,it takes some memory,but doesn't release that much, meaning on every function call a memory leak happens. **The program has crashed when calling this method about 8,000 times.**Following is the minimal code version of this behaviour:

public static bool ExportLasFromAE(ILasDatasetLayer lasdLayer, IGeometry geometry, string lasSavePath, bool overwrite)
{
    using (ESRI.ArcGIS.ADF.ComReleaser pComReleaser = new ESRI.ArcGIS.ADF.ComReleaser())
    {
        ILasDataset pLasDataset = lasdLayer.Dataset;
        ITrackCancel pCancel = null;
        ILasFilter pLasFilter = new LasFilterClass();
        pComReleaser.ManageLifetime(pLasFilter);
        pLasFilter.AreaOfInterest = geometry;
        ISpatialReference pSpatialReference = lasdLayer.Dataset.SpatialReference;
        string pFileDirectory = System.IO.Path.GetDirectoryName(lasSavePath);
        string pFileName = System.IO.Path.GetFileNameWithoutExtension(lasSavePath);
        string pFileExtensionName = System.IO.Path.GetExtension(lasSavePath);
        //IEnvelope pGeometryEnvelop = geometry.Envelope;
        IRelationalOperator pRelationalOperator = geometry as IRelationalOperator;
        for (int i = 0; i < pLasDataset.FileCount; i++)
        {
            ILasFile pLasFile = pLasDataset.get_File(i);
            IEnvelope pLasFileEnvelope = pLasFile.Extent;
            if (pRelationalOperator.Disjoint(pLasFileEnvelope))
            {
                continue;
            }
            int pFileIndex = i;
            string pFilePath = string.Format("{0}//{1}_{2}{3}", pFileDirectory, pFileName, i.ToString(), pFileExtensionName);

            if (File.Exists(pFilePath))
            {
                if (overwrite)
                {
                    File.Delete(pFilePath);
                }
                else
                {
                    continue;
                }
            }
           pLasDataset.Export(pCancel, pFileIndex, pLasFilter, pFilePath, pSpatialReference, false, false);
        }
    }
    return true;
 }


texts: Exporting climatic tabular data stored in TIF file into CSV

texts:
I have a TIF file containing tabular climatic data, described here, which I am trying to export to spreadsheet / CSV format:
http://www.helioclim.org/linke/linke_helioserve.html

In summary, it is a TIF with no header, 2160 rows by 4320 columns, where the upper left corner is 90 N, 180 W, and the lower right corner is 90S 180 E. Standard stuff. I also know that the output data is supposed to be an integer roughly between 40 and 140.

All I want to do is to extract the tabular data into a spreadsheet as source data for calculation purposes (I'm well aware that the column count exceeds Excel's capabilities). It seems the TIF is used as a very efficient container of the GIS raw data. And since the TIF contains data that is already in tabular format, I thought it would be relatively easy (not the case).

So far the only thing that has 'kinda' worked I have tried freeware "Grid Convert" to CSV, but it produces 3 columns of output data, the first 2 being all zeros and the third, 7 digits long, presumably 9,331,200 rows. It's still useless, but nothing else has worked. I don't have much experience working with this kind of data format, a starting heads up would be welcome :)

In case it might help, this is a download link to part of the raw data, under "Linke Turbidity Factor: Monthly maps"
http://www.soda-is.com/eng/services/climat_free_eng.php

Anybody have any ideas how I may extract the raw data from the TIF file?

Thanks in advance


texts: Can't seem to get QGIS 3.2 installed in Linux Mint - Only 2.18

texts: I'm at the end of my rope here.  I'm trying to install QGIS 3.2 on my computer running Linux Mint 18.3 and I've gone through installing the unstable PPA and updating and upgrading half a dozen times and it just installs 2.18.  I'm supposed to be using 3.2 for a class, and I don't understand why I'm getting 2.18 when I am not using the LTR PPAs.  If someone could give me an idiot-level explanation of how to get 3.2 installed (and more important WHAT I am doing wrong), it would be a lifesaver.


texts: Incorporation the Ordnance Survey map to Sheetsee.js plugin for Leaflet/Mapbox

texts:
I am currently using this plugin:

https://github.com/jlord/hack-spots

which works really nice in conjunction with Google Sheets.

My goal is to convert this site into my own purposes. One of the step is the change of map tile, which belongs to Mapbox in this case.

As we know, Leaflet has only 1 crs. I believe, that Mapbox too.

I tried to incorporate the Ordnance Survey map layer as per the plugin here:

https://github.com/OrdnanceSurvey/OS-Maps-API

using code such this:

       var epsg27700 = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +datum=OSGB36 +units=m +no_defs";
var crs = new L.Proj.CRS(
'EPSG:27700',
epsg27700, {
    transformation: new L.Transformation(1, 238375, -1, 1376256),
    resolutions: [896.0, 448.0, 224.0, 112.0, 56.0, 28.0, 14.0, 7.0, 3.5, 1.75, 0.875, 0.4375, 0.21875, 0.109375],
});

    document.addEventListener('DOMContentLoaded', function() {
        var gData
        var URL = "1heE0vbbWsQD3Gd2VyAYKukzyyE4sduoJ4dH5T_mWQUA"
            Tabletop.init({ key: URL, callback: showInfo, simpleSheet: true
    })
    })

    // so long, so messy
    function showInfo(gData) {
        tableOptions = {
            "data": gData,
            "tableDiv": "#hackSpotsTable",
            "filterDiv": "#tableFilter"
        }
        // make the table, and the search bar
        Sheetsee.makeTable(tableOptions)
        Sheetsee.initiateTableFilter(tableOptions)

        // create geoJSON with coordinates and other
        // useful bits from the original data
        var optionsJSON = ["order", "address", "city", "status"]
        var geoJSON = Sheetsee.createGeoJSON(gData, optionsJSON)

        // create map, tilelayer (map background), markers and popups
        var map = Sheetsee.loadMap("map", {
        crs: crs,
        layers: Road,
        zoomControl: true,
        maxZoom: 13,
        minZoom: 0,
        })

       Sheetsee.addTileLayer(' https://api2.ordnancesurvey.co.uk/mapping_api/v1/service/zxy/{tilematrixSet}/{layer}/{z}/{x}/{y}.{imgFormat}?key={key}', {
key: OS_API_KEY,
tilematrixSet: 'EPSG:27700',
layer: 'Road 27700',
imgFormat: 'png',
continuousWorld: true
     })
                   var markerLayer = Sheetsee.addMarkerLayer(geoJSON, map, "

{{ order }}", true)

        // find the latest spot with the most important
        // info filled in (to prevent map breaking if
        // someone is currently editing spreadsheet)
        var theLatestSpot = findLatestCompleteSpot(gData)
        var latestSpot = Sheetsee.ich.latestSpot({
            rows: theLatestSpot
        })
        // set it and pan to it
        $('#latestSpot').html(latestSpot)
        map.setView([theLatestSpot.lat, theLatestSpot.long], 14)

which is implicit, because the layer comes from Leaflet to Mapbox.

Unfortunately it looks like there is no option to get OS Api for Mapbox (only for Leaflet and OpenLayers)
https://astuntechnology.github.io/osgis-ol3-leaflet/ol3/01-OSGB-MAP.html

apart from this plugin, that requires Python code invoke.

https://github.com/OrdnanceSurvey/os-mapbox-adapter

Can anybody help me to solve it? How can I incorporate Ordnance Survey to this map?


texts: Process to calculate distance from points to line in QGIS

texts:
I am trying to calculate the distance between each point(as seen in the picture below) and its nearest line which is a road.
I have tried to rasterize the data on all the roads and I end up with the rasterized output which can be seen in the second image. I am not sure why the output has those numbers.


texts: Is such way to import the NetCDF file in R correct?

texts:
Recently, I asked the same question on the Stack Overflow forum and was directed to this forum: https://stackoverflow.com/questions/65467765/is-such-way-to-import-the-netcdf-file-in-r-correct.
To sum up, I have a CMORHP product (satellite based precipitation data) that I would like to import into R and extract point values for specific locations. An example, of the NetCDF file and stations is uploaded here: https://transferxl.com/08jzws5Rp8GPGp.
I am using next code to import the data into R and extract point values for these locations (stations):
library(rgdal);library(ncdf4);library(raster)
nc <- nc_open("G:/CMORPH/CMORPH_V1.0_ADJ_8km-30min_1998050101.nc",verbose=TRUE)
cshp <- readOGR(dsn="G:/CMORPH",layer="stations")
brick("G:/CMORPH/CMORPH_V1.0_ADJ_8km-30min_1998050101.nc")
data <- ncvar_get(nc, "cmorph")
data.slice <- data[, , 1]
lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
lon[lon > 180] <- lon[lon > 180] - 360
dat1 <- list(lon,lat,cbind(t(data.slice)[,2475:4948],t(data.slice)[,1:2474])) ## This could be merged with the line below
r <-raster(x=dat1[[3]],  xmn=-180, xmx=180,ymn=-60, ymx=60,crs=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
r <- flip(r, direction='y')
neki <- extract(r,cshp)

If I compare (visually)
plot(r)

and visual representation using Panoply viewer (https://www.giss.nasa.gov/tools/panoply/), both plots look similar. Moreover if I compare r and brick("G:/CMORPH/CMORPH_V1.0_ADJ_8km-30min_1998050101.nc")outputs, main properties look fine. However, if I extract data for longer period (i.e. multiple months) and compare the extracted data (precipitation) with measured point data the difference is much larger as I would expect.
The code works fine but I am not completely sure if the import of the NetCDF file into R is done correctly. The main reason for my doubts is because extracted data from multiple .nc files differs significantly from measured point precipitation. Thus, I think that I could be doing something wrong when defining projection properties. Any feedback?


texts: ArcMap - error 000728 field does not exist within table

texts: I'm trying to join the table from Excel (my calculation about one region in Poland) to my shapefile (which contains whole Poland map and different calculations, so I need to attach a new table to be able to work on my region) and run grouping analysis on it. But there's a problem with a unique ID field - it changes its name so that the grouping analysis can't be performed. It says error 000728: field does not exist within the table. Field UID, which was created by my professor, should be still named UID, but after joining my table it changes to 2005.UID. Actually all field names are getting prefix 2005. My professor has never seen anything like it so he can't help me. Why it modifies when it shouldn't? Did someone have a similar issue before?


texts: "Custom" Widget not showing in WebAppBuilder Developer

texts: I customized the in house WAB Filter widget to do a few additional things. I did not want to lose the standard filter so I just changed the folder name and kept the In House filter widget as is. For some reason, the new Filter widget I customized is not showing up. I changed the "name" as seen below. Also changed the name in Manifiest.json. Is there something else I would need to change to make it a "custom" widget to show up?


texts:
The image below shows the problem I'm facing.
I'm trying to align a DEM with the ArcMap Imagery Basemap.
Unfortunately, I can't make them to perfectly align.
I tried changing the coordinate system to WGS_1984_Web_Mercator_Auxiliary_Sphere as I understand that is what the Basemap is using but without luck.
The DEM originally uses WGS_1984.


texts: Fix the data type properties of fields when a shapefile is manipulated in R

texts:
I manipulated a shapefile in R with this code:

d<-readOGR("RECFE.shp")
dsub<-d[d$MUNICIPIO==86,]
writeOGR(obj=dsub, dsn="data", layer="ledes", driver="ESRI Shapefile")

In the dsub dataset some fields have not the same properties as RECFE.shp dataset. The precission of double fields is higher, that is, it has got more number of decimals.

       RECFE.shp                      dsub.shp
TYPE    LENGTH PRECISION   |    TYPE    LENGTH   PRECISION
DOUBLE   23        0       |   DOUBLE     23          15

Is it possible to mantain the same fields properties of the old shapefile in the created new shapefile?.


texts: How do I bring the Z value to attribute table?

texts:
I have a shapefile exported from a hydraulics software that I use. The file is exported as polylines. I convert it to a point layer using the "Extract Nodes" function. The remaining piece of information I need is the elevation (z value) of each point. I know they exist because I can see them, in the Vertex editor and the wkt_geom column when I copy paste to Excel. I see the "z()" function but can't figure out the appropriate syntax.

I tried this: z(geom_from_wkt( [shapefile_name]))

Is there a function, or a work-around, similar to the $x and $y functions to extract the z value from the geometry and display it in the attribute table?

I am talking about a point shapefile, and using QGIS 2.18.


texts: Geocoding Levels - zip code versus block face

texts:
I'm looking at some massive data, most of which has been geocoded to the block face (front of the building) level, while some of which has been geocoded to the zip code level (less accurate), giving me a lat/long pair for every sample. When we run a test on the accuracy of the zip code level geocoding, there's a largish error compared to address level coding (test was possible since some limited sample addresses were available even for samples coded at zip code level) - a median difference of 0.63 miles.  However, on plotting the lat/long pair on google maps, I see a very specific location.

What then does the error tell us?
Is the lat/long point plot (from zip geocoding) an erroneous one?
Could anyone explain the difference between zip and block face geocoding, in the sense of loss of accuracy for the former?
Finally, if I am to convert lat/long and zip code information to census tract, does the fact that the sample has been coded to zip level make a precision difference in this exercise?

Many thanks!


texts: Loading custom prj conflicts with existing CRS

texts:
I'm using ESRI Basemap (EPSG:3857) in QGIS. I downloaded some government data but the .prj-file (Clarke_1866_Transverse_Mercator) is not being used. Instead, QGIS loads in EPSG:2016 and it's not showing the layer.

When I try to create a custom CRS using the data from the .prj-file, QGIS prevents me from stating it is the same as EPSG:2016.

Here is the data .prj-file:

PROJCS["Clarke_1866_Transverse_Mercator",
GEOGCS["GCS_Clarke_1866",DATUM["D_unknown",
SPHEROID["clrk66",6378206.4,294.9786982]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["false_easting",304800.0],
PARAMETER["false_northing",0.0],
PARAMETER["central_meridian",-79.5],
PARAMETER["scale_factor",0.9999],
PARAMETER["latitude_of_origin",0.0],
UNIT["Meter",1.0]]


texts: How do I set the highlight color of a result set?

texts: I'm working with Geocortex Essentials configuring a viewer based on the Viewer for HTML 5.2.2. My workflow is functioning as designed, confirmed by the correct placement of pushpins. What I really want is the result set from my workflow to be highlighted on the screen (yellow would be great). How do I change the highlight color in the HTML5 viewer? Is this functionality the Silverlight viewer has but HTML5 viewer currently lacks?


texts:
I'm exporting from QGIS 3.4 to ACAD .dxf format.
Since I made this in QGIS 2.18 all features exported were blocks in Autocad, but now I find out that every single geometry is exported exploded in the .dxf format.

How do I export them as blocks?


texts: Seeing remote directory of files when using local ArcGIS server manager session to remote server?

texts:
Here is my setup -

AWS cluster for my ArcGIS Server instance with SD files copied to it

Currently I have to open a remote session and launch the AGS manager from the cluster so I can point to the files I copied over.

What I would like is to be able to open an AGS manager window on my local machine and point to the files already copied up. Any thoughts how I could navigate to the remote server's local paths?


texts:
I am intersecting a base road layer with 6 separate road layers to create 6 line feature outputs. I am using Model Builder for this.
I was wondering is there a way I can SUM all of the values in a field for each of the separate road layers output?
My other option is to export the road output attribute tables to excel and make the calculation there but I'd much rather if I could add such a calculation to model builder.


texts: Count how many times an attribute in a polyline intersects a polygon in QGIS

texts:
I am trying to count how many times a polyline with a specific attribute intersects a polygon in QGIS. For example, I have a polyline with a feature class called "size", this contains a text value i.e. small, medium & large. I have also made a grid of 100m x 100m squares.

I would like to be able to count how many "large" polylines intersect each grid square polygon, and then repeat this for "medium" etc.


texts: GWC tiles getting deleted for grouplayer when i change one layers style that is associated with group layer

texts:
can anyone face issue of GeoWebCache (GWC) tiles automatically getting deleted when a layers style is changed that is associated with group layer. Normally it should not happen like this., but i have tested 2 times and confirmed.

Anyone suggest me for useful documentation, where i can get complete indepth information about Geoserver gwc processing?


texts: Trying to resize a collapsed legend item in QGIS

texts:
I would like to resize an item which seems way too small in my legend (see picture below)

I would like to resize it so it looks like this:

I don't understand why it scales like this when I put it as a legend to be honest, if you have an answer to that question, I would be thankful aswell.


texts: Plotting point data from one dataset on polygon

texts:
I have a point dataset with information on: date, number of beneficiaries, and type of assistance.

I want to plot/map this info on a country polygon. I tried joining attributes by location, but that didn't work for the text fields.

Any ideas on how to map the point data on the polygon?


texts: Problem importing sqlite into Postgres (with PostGis) using ogr2ogr

texts:
I am trying to import a very large sqlite file into Postgres SQL 9.3.

The first tool I tried to use was OGR2GUI (http://www.ogr2gui.ca/). The 64 Bit version does nothing, whilst the 32 bit version says

unable to open ogr2ogr !

I have tried running org2ogr from the command line too, I have tried:

ogr2ogr --config SQLITE_LIST_ALL_TABLES YES -f "PostgreSQL" PG:"dbname='mydatabaseName' active_schema=myschema schemas=public host='localhost' port='myportnumber' user='myusername' password='mypassword' " -overwrite -lco SCHEMA=myschema -lco GEOMETRY_NAME=geometry -lco SPATIAL_INDEX=YES myfilename.sqlite  myfilename

ogr2ogr -f "PostgreSQL" PG:"dbname=mydatabaseName" myfilename.sqlite myfilename

ogr2ogr -f PostgreSQL PG:"host=localhost user=myusername password=mypassword dbname=mydatabaseName" -lco LAUNDER="YES" myfilename.sqlite -skipfailures

And various iterations on these, each time I get:

Unable to open datasource'myfilename.sqlite' with the following drivers
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> SDTS
  -> TIGER
  -> S57
  -> DGN
  -> VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> GML
  -> GPX
  -> KML
  -> GeoJSON
  -> GMT
  -> WAsP
  -> PostgreSQL
  -> PCIDSK
  -> OpenFileGDB
  -> XPlane
  -> AVCBin
  -> AVCE00
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> PGDump
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> PDS
  -> HTF
  -> AeronavFAA
 -> EDIGEO
 -> SVG
 -> Idrisi
 -> ARCGEN
 -> SEGUKOOA
 -> SEGY
 -> PDF
 -> SXF

Of course sqlite is not on this list, but how can I add it?

I have tried running osgeo set-up, but there is no option I can find to add it.


texts: MultiPointZ wkt geom "multiple part to single parts" is losing Z coordinate

texts:
I have a "MultiPointZ" unique wkt geom feature constituted of several hundreds of points (X,Y,Z) in a shapefile.

To split all these points each in a single feature, I used the "Vector > Geometry Tools > Multipart to singleparts".

It works fine except that it ended up only with "Point" geom features instead of "PointZ", thus losing the Z coordinate.

Any ideas about how to fix that?


texts: Converting ArcMap to ArcGlobe code using ArcPy

texts:
I have a code that should read the mouse coordinates, but it only works in ArcMap, can I somehow translate it into code for ArcGlobe, otherwise it gives an error: CreateObject cannot open map document.
def onMouseDownMap(self, x, y, button, shift):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = mxd.activeDataFrame
        sr = df.spatialReference
        point = arcpy.Point(X=x, Y=y)
        pointGeometry = arcpy.PointGeometry(point,sr).projecntAs(arcpy.SpatialReference(4326))
        X=pointGeometry.firstPoint.X
        Y=pointGeometry.firstPoint.Y
        pythonaddins.MessageBox("Your mouse clicked at " + str(X) + " , " + str(Y),"My Coordinates:")


texts:
One question to Archive mode in ArcGIS.
I am new to ArcObjects. Is there the possibility to change an archive class table with arcobjects by listening through events. Every time when I make an edit, the archive class is updated automatically but not with necessary information I want.  So if editing has stopped, a column to archive class table should have been added and filled with information like 'inserted', 'deleted'.

In generally I have an ArcGIS add-in button callling a userform when pressed.
So one class button.vb, one class form.vb
Where should I implement events overall. Events like mouseover, editingstopped or any other event from ArcObjects.


texts: ArcGIS javascript api application page shows blank page in Internet Explorer

texts:
I have a ASP.net application from where I send query parameters (using hyperlink) to a ArcGIS javascript application to show query result on map. Now the issue I am having is, when I first send the query parameter, the map displays correctly, but from second time it starts to show blank page. But then if I refresh IE, it shows the map correctly.

Here are some more information.

 1) I am sending different query every single time.

 2) It works perfectly in Firefox and Chrome.

 3) I don't get any error in IE console.

 4) I tried using fiddler to monitor and it looks like browser is not doing anything.

Has anyone ever encountered this issue?


texts: How to determine the columns generated when uploading a KML into a Google Fusion table

texts:
I created a Road Segment shapefile with ArcGIS desktop, the shapefile has the following columns:'CellID' 'RoadName' 'Description' Like so:

I converted into a KMZ file which I then converted into a KML. My goal is to upload that kml to a google fusion table and display it on a google map and merge it with the rest of its data based on the CellID column, when I upload the file to a Fusion Table there are only 3 Columns present Like so:

How can I get the 'CellID' column to generate in the Fusion Table? my real data has more columns however If I cannot figure out how to do it with basic data.


texts: Plotting data without geographical coordinates besides the main map in ggplot2

texts:
I am trying to plot some data as a bubble chart on the world map using ggplot2 package in R. I am able to proceed as follows.

# Load necessay packages
library(ggplot2)
library(rworldmap)
library(rgdal)

# Prepare the base world map
data(countriesLow)
world <- countriesLow
rm(countriesLow)
world <- spTransform(world, CRS("+proj=robin"))

cent <- data.frame(id = world$ISO3,
                   coordinates(world))

world <- fortify(world, region = "ISO3")

# Prepare sample dataset
data <- data.frame(id=sample(x=unique(world$id),size = 10, replace = F),
                   v1=rnorm(10, 100, 50))

data <- merge(data, cent, all.x=TRUE)

# Plot the map
theme_opts <- list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         panel.background = element_blank(),
                         plot.background = element_rect(fill="cornsilk"),
                         panel.border = element_blank(),
                         axis.line = element_blank(),
                         axis.text.x = element_blank(),
                         axis.text.y = element_blank(),
                         axis.ticks = element_blank(),
                         axis.title.x = element_blank(),
                         axis.title.y = element_blank()))

P <- ggplot() +
  geom_polygon(data=world, aes(long,lat, group=group, fill=hole), fill="black", color="cornsilk") +
  geom_point(data=data, aes(X1, X2, group=NULL, fill=NULL, size=v1), color="brown1", alpha=I(7/10)) +
  scale_size_continuous(range = c(5, 15)) +
  coord_equal() +
  theme_opts +
  scale_fill_manual(values=c("black", "white"), guide="none")

Now if I have an additional data record with unknown geographical coordinates, it is not plotted as expected.

data$id <- as.character(data$id)
data[11,1] <- "UNK"
data[11,2] <- 100
data$id <- as.factor(data$id)

P2 <- ggplot() +
  geom_polygon(data=world, aes(long,lat, group=group, fill=hole), fill="black", color="cornsilk") +
  geom_point(data=data, aes(X1, X2, group=NULL, fill=NULL, size=v1), color="brown1", alpha=I(7/10)) +
  scale_size_continuous(range = c(5, 15)) +
  coord_equal() +
  theme_opts +
  scale_fill_manual(values=c("black", "white"), guide="none")

Now, I would like to plot the data with the unknown coordinates somewhere outside the map area near the legend as a single bubble. How to achieve this using ggplot2?

What I am looking for is something like this


texts: I've installed qgis_dev 3.9.0-55 through the OSGEOW64 installer on my win7 machine. I've got the error message proj_6_2.dll is missing (qgis_dev with grass) or "could not load qgis_app.dll" (qgis-dev without grass). I've uninstalled all and reinstalled it, but always the same...any hints?


texts:
I have a geotiff image from Paris metropolitan area I have generated from MODIS data in Google Earth Engine with the following projection
data_.rio.crs
CRS.from_wkt('PROJCS["MODIS Sinusoidal",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]')

Now I have the shapefile of Paris metropolitan area. I would like to set the right crs=myCrs for the shapefile but I do not know which one it is right.
# Get images from Paris
gdf = gpd.read_file('paris.shp')     ## Paris metro area
gdf = gdf.to_crs(epsg=myCrs)


texts: How to apply the same style to different columns?

texts:
For many of my datasets I have the same imported data duplicated numerous times, each iteration displaying the data values for a different column (in my case it's fish catch weight per gridded cell per year, so lots of iterations of the same grid, with a different pattern of values per year).

For presentation purposes, it strikes me as helpful to use the same scale for each year, so I go to the original data and divide the maximum value into (arbitrarily) 10 classes and manually input the values and labels after setting the desired colour ramp and input data field:

This is extremely tedious, not just the manual excel calculation part, but also the input interface for the values & labels*, hence my intention to make use of the 'save style' option. HOWEVER: when I save the style as QML and load it for the next layer, it saves everything - including the data field (in my case, year). When I change this field, QGis automatically generates the values & labels based on the data for that field, ditching all of the values & labels I set.

Is there any way to avoid this happening & thereby duplicate class number, values and labels for each field by simply loading a style?
(somewhat meta) Am I going around this the wrong way, OR, is there a different face I could use to climb this mountain?


texts: Using OpenLayers and Georeferencer plugin in QGIS?

texts:
Let me start by saying that I am total noobie at QGIS, so this might be obvious to everyone else.

I spent hours yesterday trying to georeference a map. I needed to georeference it against vegetation features, so I thought it would be easiest to bring in Google Earth imagery via OpenLayers and then use the Georeferencer plugin. But it didn't work.

My solution:

I finally figured out that when you bring in a layer with OpenLayers, it changes the Coordinate Reference System to something else (using UTM?) so that when you use the Georeferencer panel and click in the map, it brings in numeric entries that are absurd if interpreted as ° Lat/Lon. After you load your open layer, you need to restore the reference system.

Perhaps this is clear to all you professionals, but it took me a long time to figure this out.

I would post screengrabs, but this system won't let me.


texts: Select features with same value in a given column

texts:
I first tried writing this post in generalities but realized it was more confusing that just giving the long version, so here goes.

I have a tree inventory that has property quads (in front of a house, to the left, etc) and then a sequence number so that each tree can be identified by its relative location on that property. Also I have street addresses and names among other attributes. My plan is to automate a way to update the sequence numbers.

My current plan to is concatenate a field for each unique location ("Address" + "StreetName" + "Quad".) I will also create another field with those but also the "Sequence_Num" added onto it.

My plan then is to Sort but the second unique field and then create a for loop that selects all trees by the first unique field and then generates an autoIncrement() value.

The problem is that right now I might have several trees with the same exact value for the second unique field (same "Sequence_Num") and I can't be sure it will sort correctly. So I want to find all trees that have the same Unique value for the 4 field concatentation (address, streetname, quad, sequence.)

Short recap, is there any way to find all features that have the same value for a given field?

For ArcGIS.


texts:
I'm trying to export 150 individual TIFF files from a large dataset. I got a script that it seems to work but it generates TIFFs at default resolution but not at the desired one.

First, some words about my problem.
Let’s say I have a big dataset of building blocks of a city and I want to perform a batch georeferenced TIFF export based on a grid that I have constructed.
My grid is called Grid_file.shp and I wrote a script that zooms to every part of the grid and then performs the export process. It seems to work but the generated TIFFs are at the default resolution (96dpi) and not at the desired one (300dpi) any ideas?
Here’s my code:

mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, "Layers") [0]
lyr = arcpy.mapping.ListLayers(mxd, "Grid_file", df)[0]
for x in range (0,12):
 arcpy.SelectLayerByAttribute_management(lyr, "NEW_SELECTION", ' "Id" = %d '%(x))
 df.zoomToSelectedFeatures()
 arcpy.mapping.ExportToTIFF (mxd, r'C:\temp\%d.tif' %(x), df, 5030, 2855, 300, True, "8-BIT_PALETTE", "NONE",True)

I'm using ArcInfo 10.2.2, any ideas?


texts:
I've spending last two weeks trying to discover how to do a Manhattan Voronoi on QGIS.
For the difference between Manhattan Voronoi and "normal" Euclidean Voronoi, please check this link

Can someone point me to the right direction?  In the link bellow you can find a sample of the information I've got (points and one polygon)

https://drive.google.com/file/d/0B4WB0ixNBUYkREtlcW0weDJUYU0/edit?usp=sharing

I would like to do in on QGIS, but if it's not possible, I'm open to other solutions.


texts: QGIS - How to render QgsMapCanvas on QApplication's main window

texts:
I am working with C++, QT 5.11.3, Visual Studio 2015 and QGIS 3.7

I'm trying to develop QT c++ application that uses the qgis api. I want to display simple raster and vector map in my qt application.
I create QgsVectorLayer and QgsRasterLayer that added to the QgsProject instance and then QgsMapCanvas object. That is no problem so far.

But when i turn the visible property of QgsMapCanvas object, the application crash. The error message and the crash window is above.

  QCoreApplication::sendEvent: "Cannot send events to objects owned by a different thread. Current thread ... Reciever (of type 'QGSMapCanvas') ...

I try to this process at different qt thread or main thread no matter what, does not help.
What can i do wrong?


texts:
I am trying to apply a code (in R) for all images that I have in a folder.
The code works for 1 image, but when I try to make a loop out of it, it gives me error.

library(raster)
library(rgdal)
ET<-"ET"
pre<-"pre"
files <- list.files(path="......", pattern="*.tif", full.names=T, recursive=FALSE)
for (r in files) {
name = r
subname<-sub(ET,pre,r)
srraster<-stack(r,subname)

#identifying some pixels as NA
r[r>30000] <- NA
# creating a function for replacing NA pixels
fill.na <- function(t) { center = 0.5 + (width*width/2)
  if( is.na(t)[center] ) { return( round(mean(t, na.rm=TRUE),0) )  }
else {return( round(t[center],0) )}}
width = 15
r= raster(focal(r, w = matrix(1,width,width), fun = fill.na, pad = TRUE, na.rm = FALSE))
#Saving files
savename<-sub("res",ET,basename(file))
  writeRaster(r,file=savename)

}


texts:
I have created a simple mapboxgl map, and I am adding a source that contains point data. I am trying to then apply a style to the data. The approach I am taking works well for polygon & line data, but I can't find the magic cartocss that will render points. Here is the example: http://jsfiddle.net/4q1db482/

The reason I do not want to simply use the getJSONSource method for point data is to avoid bringing a large number of vectors into the browser.

The question is two part:
1. How can I apply a style to this data that I have created in Mapbox-studio?
2. Or how can I create a style that will work for some simple point data and push it into the style.layers array? It's a bit cumbersome that 'marker-line-color' cartocss is valid in mapbox studio, but not in mapboxgl. I am trying this:

   style.layers.push({
        'id': 'dataSource',
        'source': 'tom',
        'interactivity': 'true',
        'source-layer': 'us-ski',
        'render':{
            'text-field': '{name}',
        },
        'style': {
            'text-color': '#8786e0',
            'text-font': 'Open Sans Semibold, Arial Unicode MS Bold',
            /** This is the style used in mapbox-studio **/
            //"marker-symbol": "monument",
            //'marker-line-color': '#3333FF',
            //'marker-width': '3',
            //'marker-fill': '#8786e0'
        },
        'type': 'symbol'
    });


texts: Giving input of multiple rasters in Raster Calculator of ArcGIS Desktop

texts:
I want to add-up more than 100 raster datasets to derive the sum of all.
Adding all these one by one in calculator using basic algebra expression is time consuming.
Can anybody suggest me how can I add-up all these rasters using a short expression?
You can suggest me any other tool in ArcGIS Desktop too.


texts: How to set the extent of a PostGIS geography table?

texts:
I have multiple geography tables in PostGIS database. All of them have extents of (-180,-90;180,90), probably in their metadata, which makes editing difficult. I'd like to set the extent to exactly the bounding box of them, but I can't find any command to alter them.

Extents are viewed in QGIS by right clicking the layer and click Properties, then click Info tab.

How do I set the extents to exactly the bounding boxes of my geography tables?


texts: VETiledLayer in ArcGIS JavaScript API v4?

texts:
I'm a first-time arcgis js user, and I'm trying to figure out how to use bing maps (or google maps) with it. So far I found out that with v3 you used VETiledLayer – however that seems to be gone in v4.

what's the way to do this in the latest version?


texts:
Could somebody tell me what is wrong with this code? This is a part of my script and returns totally different result from "Nearest Hub" with the same line and point layers. I would like to identify the nearest link to each node.

lineLayer = QgsMapLayerRegistry.instance().mapLayersByName("road links_link")[0]
            spIndex = QgsSpatialIndex() #create spatial index object of links
            feat = QgsFeature()
            feat1 = QgsFeature()
            fit = lineLayer.getFeatures()

            # insert links to index
            while fit.nextFeature(feat1):
                spIndex.insertFeature(feat1)

            while stp.nextFeature(feat):
                pt=feat.geometry().asPoint()
                nearestIds = spIndex.nearestNeighbor(pt,1) # from stop to link
                featureId = nearestIds[0]
                #gets the attributes of the link
                fit2 = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId)).next()


texts: Visualize NetCDF4 file containing burnt area data with QGIS

texts:
I've been trying for a while to get QGIS to visualize a ~300 MB .nc file (netCDF4).
Data is burnt area data from here: Copernicus Global Land Service.
What I currently do is either drag the nc file to QGIS, after which it imports the 4 layers with the correct names, or Layer > Add layer > Raster layer, which results in the same correct import.
However, I can't find out how I should get QGIS to visualize the data for me. From what I can tell from online sources, the layer should be shown right after you add it. For me it shows this:
There is no error or anything, it just doesn't seem to visualize my data file. In the bottom right it says Unknown CRS, clicking that and setting WGS84 doesn't change it from saying Unknown CRS.
Does anyone know what I'm doing wrong here? Does the file just have nothing to show? I selected south Portugal in August 2018, when viewing in "quick look" on the Copernicus site there is some burnt area data there.


texts: Is there a standard vertical datum for meters below the WGS84 ellipsoid, or below MSL?

texts:
I have data I am publishing a WMS layer representing data under the surface of the water, say temperature at 100m depth.  What is the accepted way to represent this in WMS GetCapabilities document?

Do I define the dimension like:

<Dimension name="elevation"
           units="EPSG:5030"
           unitSymbol="m">0,50,100</Dimension>

or

<Dimension name="elevation"
           units="EPSG:5030"
           unitSymbol="m">0,-50,-100</Dimension>

or something else or should a custom dimension be used in this case?

EPSG:5030 is a length parameter, it doesn't indicate that it's related to any projection or that it might indicate distance from mean sea level for that matter. see http://trac.osgeo.org/geotiff/wiki/VerticalCS
EPSG:5030 represents height above the WGS85 ellipsoid. I found a reference to EPSG:5715, but it doesn't say much other than it's 'msl depth' - does this mean values below msl should be positive?


texts: ArcPy FieldMappings syntax using Merge tool

texts:
I'm working on script that takes several table results from the Zonal Statistics As Table tool and merges them into a single table. When merging I encounter the issue that the exported .dbf files have different lengths for a field I need, "UID". The merge takes the first table added and uses this length as the max for the merged "UID" field.

What would the syntax look like to achieve this?

Here is my code:

arcpy.Merge_management(listTable5M,outTable+"\\"+"5M_zonal_statistics_elevation_merged.dbf","UID \"UID\" true true false 25 Text 0 0,first,#")

Do I need to use the FieldMappings() object?

The first answer makes a lot of sense. I'm still not getting the result I want however.

The Merge still throws an error. Here is your code implemented into the script:

dbf_list_5m = arcpy.ListTables()
    arcpy.AddMessage(dbf_list_5m)
    max_UID_length = 0
    field_mappings = arcpy.FieldMappings()
    for table in dbf_list_5m:
        fields = arcpy.ListFields(table)
        arcpy.AddMessage(table)
        for field in fields:
            if field.name == "UID":
                if field.length > max_UID_length:
                    max_UID_length = field.length
                    arcpy.AddMessage(max_UID_length)
                    field_mappings.removeAll()
                    field_map = arcpy.FieldMap()
                    field_map.addInputField(table, field.name)
                    field_mappings.addFieldMap(field_map)
    arcpy.Merge_management(dbf_list_5m,outTable1+"\\"+"5M_zonal_statistics_elevation_merged.dbf", field_mappings)

It is getting the max_UID_length (17 in this case) but still bombs out during the merge. I can successfully merge using the gp tool if I change the length manually (or add the table with the 'UID' field 17 characters in length first).


texts: Selecting features using two expressions with PyQGIS

texts:
I am selecting a river by expression. There are more than two rivers with the same name. So, I am thinking of using HUC. I found the following site about the expression:
Selecting features using an expression with PyQGIS
But I could not get a valid syntax if there are two expressions. My code is the following. How do I fix it?
layer.selectByExpression("\"PNAME\"='{}'".format(river_name)AND "\"HUC\" = '{}'".format(huc))


texts: Layer legend symbol in OpenLayers from ArcGIS REST API

texts: I am using OpenLayers to display vector layers coming from Esri ArcGIS REST API. Each layer has an attribute called renderer which has the information needed for drawing a legend. Is there a built-in way or 3rd party library to parse that information to get a layer symbol?


texts: Removing all layers from group using PyQGIS?

texts:
I want to remove all layers from a distinct group.

root = QgsProject.instance().layerTreeRoot()
the_group = root.findGroup("my_group")

parentGroup = the_group.parent()

I think I have to loop the elements of the group and remove them

for i in parentGroup.children():
...

And then calling the removeMapLayer() method within the loop.

But I can't put these things together.


texts:
I have some archaeological scanned raster that they have a north arrow and a scale bar. There is no grid, no coordinate system.I would like to put those files in a general map. I don't care about the exact position but I would like to keep the scale and the north direction.  So I need some tools like move raster, rotate raster, scale raster
 Are there such tools in QGIS?
Thanks in advance for your help.


texts:
I need to import data on WAsP software but I am finding problems to convert the data I got from an official website to a file readable by WAsP.

in summary I have a file like this:

 -39.0250     -6.97502      314.000
 -39.0242     -6.97502      310.000
 -39.0233     -6.97502      316.000
 -39.0225     -6.97502      312.000
 -39.0217     -6.97502      306.000
 -39.0208     -6.97502      309.000
 -39.0200     -6.97502      318.000
 -39.0191     -6.97502      316.000

and so on.....very looong file.

Do you have any idea how to turn a .txt like this to a .xyz

I am beginner on such topic, so please, be as clear and specific as possible.


texts: Geoserver 2.13+ GetCapabilities is Missing XML Encoding

texts:
Hoping someone can help me on this.

So I've installed Geoserver 2.11 and am able to produce proper XMLs from my GetCapabilities.
For example:

<?xml version="1.0" encoding="UTF-8"?>
<Capabilities version="1.0.0" xsi:schemaLocation="http://www.opengis.net/wmts/1.0 http://schemas.opengis.net/wmts/1.0/wmtsGetCapabilities_response.xsd" xmlns:gml="http://www.opengis.net/gml" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns="http://www.opengis.net/wmts/1.0">

This works fine for me inside my ArcGIS runtime.

But when I use version 2.13+ it is missing the XML encoding at the header of the GetCapabilities which is no longer recognized by my ArcGIS runtime.
I have tweaked every setting I could think of but cannot find a solution to this.
This is what my GetCapabilities looks like in the newer versions of 2.13+.
Notice there is no XML encoding, so my ArcGIS runtime does not recognize it as a valid XML. Does anyone have any suggestions? I have lost hours on this seemingly silly thing.

<Capabilities xsi:schemaLocation="http://www.opengis.net/wmts/1.0 http://schemas.opengis.net/wmts/1.0/wmtsGetCapabilities_response.xsd" version="1.0.0"><ows:ServiceIdentification><ows:Title>GeoServer Web Map Tile Service</ows:Title>


texts: Comparison of % change in two photos with GRASS

texts:
I have two aerial photos of the same location taken at two different dates. Is there a way to compare the % difference between those two images using GRASS or is there a better tool?

Sorry, also if there are other tips and tricks I am relatively new at this.


texts: Converting from Eastings and Northings on British Ordnance Survey to WSG84

texts:
I am a software developer currently working on a project that relies heavily upon knowing if our position is within a defined region - specified as a polygon using a set of points to define a linear ring.

The points are defined using Eastings and Northings from the British Ordnance Survey OSGB36. We use a nuget (software) package called GeoUK to convert this into a WSG84 latitude and longitude using the code below:

public LatitudeLongitude returnLongLat(float flUKx, float flUKy) {

        var cartesian = GeoUK.Convert.ToCartesian(new Airy1830(),
            new BritishNationalGrid(),
            new EastingNorthing(
                flUKx,
                flUKy));

        var wgsCartesian = Transform.Osgb36ToEtrs89(cartesian);
        var longlat = GeoUK.Convert.ToLatitudeLongitude(new Wgs84(), wgsCartesian);
        return longlat;
    }

When we look at the results of this on Google Maps / Earth however, we are noticing errors in the order of up to 5 meters - usually things are a couple of metres off...

I've seen suggestions that this is a just a fact of life with this particular transformation, but that sounds like something of a cop out.


texts: Displaying label text when polygon layer is rotated in QGIS?

texts:
How can i have the text labels display when the layers are rotated? I have several polygons layers rotated and i can not get the label text to display. When the rotation is removed the label texts shows but when it is rotated again it goes off.

I am running QGIS 2.8.2 Wien.


texts:
What I'm attempting to do is create a list of .dbf tables using a wildcard which I've managed to do. From there my goal is to search the list and extract the value of a specific field. This field only contains 1 row, the value will then be used to populate a field contained in another table which, I haven't begun to attempt yet.Below is a code snippet:

tableList05 = arcpy.ListTables("*_2005.dbf")
#for table05 in tableList05:
with arcpy.SearchCursor(tableList05, ['SUM'])as cursorFive:
    for row in cursorFive:
       print row

For now I'm simply trying to get those values to print before trying to manage how to copy them and update the records in an existing table. The error message I receive when I run this section of code tells me the .dbf tables don't exist.


texts:
I am installing QGIS and GRASS as part of an FOSS GIS class I'll be teaching this summer. I used the OsGeo4W installer. This installer has always worked flawlessly for me except on the machines at the college.

Initially the install gave me a "The ordinal 3109 could not located in the dynamic link library LIBEAY32.dll" error. I remembered the work around for this from last spring. I copied the libeay32.dll from c:\osgeo4w\bin to c:\osgeo4w\qgis\apps\bin and was able to complete the install and have QGIS operate successfully.

However, when I enable the GRASS plugin, open a mapset and then try to add data I get the following and QGIS crashes.

Cannot get projection
Cannot run module
Command: c:/osgeo4w/apps/QGIS/./grass/modules/QGIS.g.info info=project

The layer itself has a defined coordinate system and the QGIS project is set
to the same coordinate system.

If I create a new Mapset it fails to import an existing dataset, either from
a loaded dataset, or via OGR. I get this error:

Cannot start module v.in.ogr

Command: c:/osgeo4w/apps/grass/GRASS--1.2/bin/v.in.ogr.exe --interface-description

GRASS itself will not open via either GUI or commandline.

I've gone into the OsGeo4W installer/Advanced and reinstalled both QGIS and GRASS without any change.

I need to get this up and running by the end of the week. Any suggestions will be greatly appreciated. All I can think of is that perhaps there are restrictive permissions on some folders. The computer is a Windows 7 64 bit machine.

Thanks!
Kurt


texts: WFS GetFeature: Is DWithin selection equal to an Intersect selection using a circle?

texts:
I can't find any WFS implementations that support circle geometries for selection. There is the DWithin spatial selection but I can't find the specification of what is actually does.

Does it select objects where the geometry is entirely within the specified distance?
Or
Does it select objects where some part of the geometry is within the specified distance?


texts: RadMap for Silverlight reading in Shapefiles

texts:
I have created shapefiles now which is basically a couple of polygons representing rooms in a building. They appear fine in QGIS, but when I try and load it into the application (RadMap in Silverlight) the screen just goes blank.

Anybody have any experience with this?

EDIT: I'm using QGIS and have tried MapWindow too and the shapefile CRS is set to the recommended 'EPSG:4326 - WGS 84'. When I drew the polygon I made sure the scale was set to 1:1 in the hope of eliminating the possibility that the shape is just to small to see when it renders in RadMap/Silverlight.

Thanks,
Jacques


texts: Converting LAS Dataset into a point feature class with x-y-z in the attribute table?

texts:
How can I convert a LAS data set (displayed as LAS point elevation in ArcGIS 10.5) into a point feature class with x-y-z information attached to it (or stored in the  attribute table)?

ArcGIS has the LAS to Multipoint tool (and then the Multipart to Singlepart tool), but the created output doesn't include any coordinates nor height information, as far as I can see.


texts: How to fix link between layer and shapefile in ArcGIS 10?

texts:
Does anybody know how to fix a link between shapefile and layer file in ArcGIS 10?

Every time I add layer file, there is a small exclamation mark next to icon and I have to define a source (path to a relevant shapefile).

Therefore, I was wondering if there is a way for ArcGIS to remember that, so I don't have to do it every time I add that layer file?


texts: For a ESRI Geodatabase Feature Class, which table holds the spatial reference?

texts:
For a shapefile, the spatial reference / coordinate system information is held in the .prj file.  I would like to know which file or table in a geodatabase holds the equivalent information for a feature class.  For example, if I have a geodatabase called "TownOfCoventry.gdb" and I use the Windows OS to drill down into the gdb folder, there are many many files for a single feature class.

a0000001.gdbindexes
a0000001.gdbtable
a0000001.gdbtablx
etc.
etc.
etc.

I would just like to be able to tell students that while a .prj file contains the spatial reference for a shapefile, a ?????? table file contains the spatial reference for a feature class within a geodatabase.


texts: Export a very large attribute table (1,5 M rows) into Excel spreadsheet QGIS

texts:
I have a very large shapefile containing rainfall data for South Africa. Each point is a single rainfall that lasted for 1 hour. I have 1,5 million points. For further analysis, I need to export it into Excel.

When I try to save the shapefile using "Save feature as" either into an xlsx file or a CSV file, the first 1550 rows are omitted, the file starts at row 1551, and the last 200'000 rows are omitted as well. Any insight on what I may be doing wrong?


texts:
I am wanting to automate the preparation of raster layers from Sentinel-2 to be used in an unsupervised classification. So far, I have been able to write the following script, which will successfully prepare the images in such a way that each band will cover the desired extent even when the extent is covered by multiple satellite image tiles:

bandset = {
    "band02": Band02,
    "band03": Band03,
    "band04": Band04,
    "band05": Band05,
    "band06": Band06,
    "band07": Band07,
    "band08": Band08,
    "band8A": Band8A,
    "band11": Band11,
    "band12": Band12
}

for bandidx, bandimages in bandset.items():

    # Generate output file
    outputFile = os.path.join(Output_folder, bandidx+".tif")

    # Merge band images together
    algorithmOutput = processing.runalg("gdalogr:merge", Band_images, False, False, 5, None)
    ret_merged = algorithmOutput['OUTPUT']

    # Reproject the layer to ensure that all layers will have consistent projection and resolution
    algorithmOutput = processing.runalg("gdalogr:warpreproject", ret_merged, "", Target_CRS, "", Target_resolution, 0, 5, 4, 75, 6, 1, False, 0, False, "", None)
    ret_reprojected = algorithmOutput['OUTPUT']

    # Produce buffer around region of interest
    algorithmOutput=processing.runalg("qgis:fixeddistancebuffer", Region_of_Interest, buffer_size_m, 5, False, None, progress=None)
    ret_buffered=algorithmOutput['OUTPUT']

    # Clip the merged band images to the extent of the desired area
    algorithmOutput=processing.runalg("gdalogr:cliprasterbymasklayer", ret_reprojected, ret_buffered, '0.0',  False, True, True, 5, 4, 75, 6, 1, False, 0, False, "", outputFile)

Then I tried performing the unsupervised classification, in the Semi-Automatic Classification Plugin in QGIS 3.0. I created a band set with the raster layers that are generated by this script and tried running a clustering of the band set using ISODATA when I received the error message "Error [46] : Error reading raster. Possibly bands are not aligned". I thought this was strange because accessing the metadata using the code block pasted below showed that all raster layers have identical size, projection, origin, and pixel size.

for bandidx, bandfiles in bandset.items():
    filename = os.path.abspath(bandfiles[0])
    dataset = gdal.Open(filename, gdal.GA_ReadOnly)

    geotransform = dataset.GetGeoTransform()
    projection = dataset.GetProjection()
    srs = osr.SpatialReference(wkt=projection)

    print("Driver: %s, %s" % (dataset.GetDriver().ShortName, dataset.GetDriver().LongName))
    print("Size: %d x %d x %d" % (dataset.RasterXSize, dataset.RasterYSize, dataset.RasterCount))
    print("Projection: %s" % (srs.GetAttrValue('projcs')))
    print("Origin: (%d, %d)" % (geotransform[0], geotransform[3]))
    print("Pixel Size: %d x %d" % (geotransform[1], geotransform[5]))

However, when I tried aligning the raster layers using the "Raster > Align Rasters" tool before making the band set and then performing the unsupervised classification, there was no error. Thus, aligning the rasters solved my problem completely!

So my question is this: How do I call the Align Rasters tool from a python script, or how would I write a script that accomplishes the same task as the Align Rasters tool?


texts: Selecting few biggest polygons within class

texts:
I have polygon shapefile with objects that have 'class' and 'area' values in attribute table. I want to select five biggest polygons in each class. I can't do it manually because of large number of shapefiles and objects in each shapefile. Is it possible with ArcGIS that I use?
If not maybe some other software can do the trick...
Thanks in advance for help.

edit:
I can use file/personal geodatabase or SQL Server geodatabase as well if this is nessesary to make that selection.
I will not be able to create python script, I would rather use some SQL selection.


texts: Drawing a directed line based on dataframe indexes?

texts:
I trying to draw a directed graph based on a geopandas dataframe with the Cartoframe library.

The Geopandas dataframe has the following structure:

  deviceId  routeLoad  longitude   latitude                  geometry
0    depot          0   20.618602  30.959910  POINT (20.61860 30.9599)
1   pointA          8   20.057890  39.516000  POINT (20.05789 30.4160)
2   pointB         11   20.047890  39.516000  POINT (20.05789 30.5160)
3   pointC         14   20.043920  39.335000  POINT (20.04392 30.3350)

Right now, I have the possibility to display these points on the map with Map(Layer(test)) but I do not have any idea how to create directed graph with CartoFrame and if it possible to do it without too much work.

The ultimate objective is to have a path going through the different observations depicted inside the geopandas dataframe.

For instance:

depot --> pointA --> pointB --> pointC --> depot

For simplicity, this is the sample code.

import pandas as pd
from geopandas import GeoDataFrame, points_from_xy
from cartoframes.viz import Layer, Map, color_category_style, popup_element
test = pd.DataFrame([{'deviceId': 'depot', 'routeLoad': 0, 'longitude': 20.618602, 'latitude': 30.959910},
                    {'deviceId': 'pointA', 'routeLoad': 8, 'longitude': 20.057890, 'latitude': 39.416000},
                    {'deviceId': 'pointB', 'routeLoad': 11, 'longitude': 20.047890, 'latitude': 39.516000},
                    {'deviceId': 'pointC', 'routeLoad': 14, 'longitude': 20.043920, 'latitude': 39.335000}])
test = GeoDataFrame(test, geometry=points_from_xy(test.longitude, test.latitude))
Map(Layer(test, popup_hover=[
        popup_element('deviceId'),
        popup_element('routeLoad')

    ]))


texts: Why can't I perform supervised classification using Landsat 8 SR images in Google Earth Engine?

texts:
This is the link to the code.
https://code.earthengine.google.com/ba342d996547c4fa8151bc99052036cf

This error message appears

var ColecaoL8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var diltroDataL8 = ColecaoL8.filterDate('2016-01-01', '2016-12-31');
var ComposicaoL8= ee.Algorithms.Landsat.simpleComposite({
collection: diltroDataL8, asFloat: true});
Map.addLayer(ComposicaoL8 , {'bands': ['B4', 'B3', 'B2'], min: 0, max:0.3}, 'Composição Imagens Landsat 8');

var classes = urban.merge(vegetation).merge(water);

var bandas = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'];

var treinamento = ComposicaoL8.select(bandas).sampleRegions({
  collection: classes,
  properties: ['landcover'],
  scale: 30});

var classificar = ee.Classifier.cart().train({
  features: treinamento,
  classProperty: 'landcover',
  inputProperties: bandas});
//print(classificar.explain());

var classificado = ComposicaoL8.select(bandas).classify(classificar);
Map.addLayer(classificado, {min: 0, max: 2, palette: ['red', 'green', 'blue']},'Classificação');

//////////////////// VEGETAÇÃO ////////////////////
var vegetacao = classificado.select('classification').eq(1);//vegetação tem valor 1
//Map.addLayer(vegetacao);

var areaVegetacao = vegetacao.multiply(ee.Image.pixelArea()).divide(1e6);

var statVegetacao = areaVegetacao.reduceRegion ({
  reducer: ee.Reducer.sum(),
  geometry: geometry,
  scale: 30,
  maxPixels: 1e9});
print ('Área de Vegetação (km²)', statVegetacao);


texts: Parse geometry from MYSQL on the client side

texts:
I am working with the existing production Mysql database, that contains some geometry (wkb, I think) as a blob.
While documentation offers to use ST_AsText(boundaries.geometry) AS geom, Servers says "FUNCTION db.ST_AsText does not exist" - I guess, MySQL is too old to have spatial transformations.

I am reading data from python, and wonder, How can I parse those bytes into geometry on the client side?


texts: How to create unit and quantity "amount per area unit" in JSR-363?

texts:
I am trying to integrate units and quantity with JSR-363.
I know and understand how to create  SI and NONSI units and convert e.g. Areas units  from one to another.

But, in a kind of geographic cartographic way, I want to create another unit: from absolute population number to a relative on e.g. population per km²

First question: In what unit can I put numbers with plane quantities numbers
like population number or amount of eggs.

Second question:
How can I create a new unit with combination with Area and amount number
without any java.lang.ClassCastException

I am using maven packages and Java 8

</dependency>
  <dependency>
      <groupId>javax.measure</groupId>
      <artifactId>unit-api</artifactId>
      <version>1.0</version>
  </dependency>
  <dependency>
      <groupId>tec.units</groupId>
      <artifactId>unit-ri</artifactId>
      <version>1.0.3</version>
  </dependency>
<dependency>
  <groupId>si.uom</groupId>
  <artifactId>si-units-java8</artifactId>
  <version>0.9</version>
</dependency>

and the example

 Quantity<Area> area = Quantities.getQuantity(4000.0,SI.SQUARE_METRE);
 Quantity<*????*> population = Quantities.getQuantity(60, *????*);

 Quantity< ???> result= population.divide(area).getValue();

The goal is to use any number for population and any number of area, where it is also possible to switch between m² or km² or any other area unit

What I also don't really understand is, how to get units and values from a calculation in a new Quantity class ala result?

==== UPDATE

with the code example and BaseUnit I can create the result.

but with my maven packages I get an Exception

 Exception in thread "main" java.lang.IllegalArgumentException: Cannot format given Object as a Unit
at tec.units.ri.format.SimpleUnitFormat$DefaultFormat.format(SimpleUnitFormat.java:692)
at tec.units.ri.format.QuantityFormat$Standard.format(QuantityFormat.java:188)
at tec.units.ri.format.QuantityFormat.format(QuantityFormat.java:125)
at tec.units.ri.format.QuantityFormat.format(QuantityFormat.java:143)
at tec.units.ri.quantity.DoubleQuantity.toString(DoubleQuantity.java:121)
at java.lang.String.valueOf(String.java:2994)
at java.io.PrintStream.println(PrintStream.java:821)

it is working, if I use the getValue or getUnit method.
But this way is less

import si.uom.SI;
import tec.units.ri.quantity.Quantities;
import tec.units.ri.unit.BaseUnit;

import javax.measure.Quantity;
import javax.measure.quantity.Area;
import javax.measure.quantity.Dimensionless;

public class App{
    public static void main( String[] args ){

    BaseUnit<Dimensionless> PEOPLE = new BaseUnit<>("people");

    Quantity<Area> area = Quantities.getQuantity(4000.0, SI.SQUARE_METRE);

    Quantity<Dimensionless> pop = Quantities.getQuantity(1000, PEOPLE);

    System.out.println(area); //not working with exception
    System.out.println(pop); //not working with exception
    System.out.println(pop.divide(area));//not working with exception

    System.out.println(pop.divide(area).getValue());
    System.out.println(pop.divide(area).getUnit());

  }
}


texts: Deleting points that fall outside a polygon

texts:
I have a polygon that is my region of interest.  I also have a database of observation points that I have imported into QGIS as a points layer.  I'd like to exclude all the observation points that fall outside my region of interest (the polygon).  Also happy to delete those within the polygon if that is easier (I have polygons of both region of interest and region not of interest if that makes sense).  Ultimately I'm just wanting to identify those within the region and export that information (at least for now).

I can do this manually by deleting each point, but is there an automated way I can do this in QGIS?

I'm afraid I haven't managed to find anything to really try, what I've found relates to clipping rasters or polygons, not points.  Any suggestions much appreciated - thanks


texts: Is there a way to listen for the project save event in ArcGIS Pro? I need to check for something before allowing or denying the save event.


texts: Project data built with incorrect spatial reference - how to correctly re-reference all data

texts:
Just discovered an interesting situation when coming to amend some MXDs and generate area calculations.
A previous consultant - when creating the first MXD and project data, set the correct spatial reference in the data frame but must have built all data creation (digitisation) off an image that was not correctly referenced.

All the data has now been created at around -104.09°, 40,118°. This was discovered when georeferencing an engineering drawing for the site and then calculating areas from features digitised and realising the scale didn't match up.

Is there an easy way to convert all data coordinates to the correct location? Re-referencing data based on a correctly referenced set of imagery would be too time-consuming for the project needs, and it's not a situation where assigning a new spatial reference is appropriate, as it has the correct spatial reference on the dataframe and the feature classes - it's just that the coordinates are out entirely.


texts: How to create a time slider using arcgis javascript that only shows a year at a time?

texts:
I am looking at this example on how to create a time slider using the arcgis javascript api. Instead of showing the data cumulative from the start date. I want to just show the data for that particular year. I have tried changing the ThumbCount to 1, but all that does is set the start time to the start time. Other than that, I can't seem to find any more information on how to get the behavior that I want.

I understand that I can use 2 thumbs on the time slider in order to get desired result. However, I am looking to see if it is possible by using only 1 thumb slider.


texts:
I am using osgeo4w shell to execute gdal_translate command for creating geospatial pdf. My command is:
"gdal_translate -of PDF -tr 0.0001 0.0001 -a_srs EPSG:3765 Granice.tif Sisacko_moslavacka_zupanija.pdf -co OGR_DATASOURCE=Extent.vrt".
The problem I am facing is in encoding. In my Extent.vrt file I have provided connection to MSSQL, and I want to show labels of Counties names. Layer in my vrt file that shows labels:

Names in my db look fine, but written in my pdf they look something like this "VorkapiÄ‡", etc.
Where can I define GDAL encoding on how will it read and write my shapefiles from vrt? I have tried adding
--config SHAPE_ENCODING "UTF-8" in my gdal_translate command but it doesn't solve it.


texts: ERROR pj_get_spheroid_defn EXTENSION POSTGIS2_96 POSTGRES-9.6

texts:
I recently had this error, and surely other people also are going to have when they update the POSTGIS and POSTGRESQL, so I want to share the solution.

OS: CENTOS 6 (x86_64)

Installed: postgresql-9.6
from this repository

yum install postgresql96-server postgresql96

installer: postgis2_96

yum install postgis2_96

problem/error when start to enable POSTGIS using psql or pgAdminIII :

psql (9.6.3)
postgres=# \connect SANTANDER
    You are now connected to database "SANTANDER" as user "postgres".
    SANTANDER=# CREATE EXTENSION postgis;
    ERROR:  could not load library "/usr/pgsql-9.6/lib/postgis-2.3.so": /usr/pgsql-9.6/lib/postgis-2.3.so: undefined symbol: pj_get_spheroid_defn

I checked the PROJ version:

$ yum list installed | grep proj
proj.x86_64        4.7.0-1.el6                           @epel

the problem is solved upgrading the PROJ

$ yum upgrade proj
Updated:
  proj.x86_64 0:4.8.0-2.rhel6


texts:
I would like to import my Garmin Topo Maps to QGIS. I want to use the Topo Map as a Layer to create Shapefiles.

Can anybody tell me if it is possible to import the Garmin maps including waypoints as a vector layer in QGIS?


texts: Reproject vector layer in oracle or PostGIS without export!

texts: I have some vector layers in Oracle and PostGIS, point, line and polygon, and i want to reproject those layer in a different projection but i don't want to export the projected layer, i want to do it on the original layer because i have related tables with the feature ID. I know about on the fly projection but i want to permanent transform the data.


texts: Editing Layout-Items after loading a QPT-Template in QGIS 3 with Python

texts:
I am facing a problem with the new Layout API while implementing a plugin for QGIS 3.
Starting point is a layout template in qpt-Format, which contains different elements like (attribute) tables, map, labels and scale-bar. This template is loaded into the Layout Manager with the code below.

During import I want to make dynamic changes to the layout elements like filling tables with new rows or editing the label text.

Any idea on how to address elements of the layout during or after the import, in order to edit the content of labels and tables?

project = QgsProject.instance()
composition = QgsPrintLayout(project)
document = QDomDocument()

# read template content
template_file = open('C:/location/of/template.qpt')
template_content = template_file.read()
template_file.close()
document.setContent(template_content)

# load layout from template and add to Layout Manager
composition.loadFromTemplate(document, QgsReadWriteContext())
project.layoutManager().addLayout(composition)

# maybe using QgsLayoutItem???
# layoutItem = composition.itemById("1")


texts: "could not get any sources: no image returned from source WMS" error in mapproxy demo

texts:
I am trying to connect my mapserver to MapProxy. First part of the map file is shown below

MAP
FONTSET "fonts.lst"
IMAGETYPE png
MAXSIZE 4000
SIZE 800 800
EXTENT -20000000 -20000000 20000000 20000000
UNITS meters
IMAGECOLOR 153 179 204
WEB
   METADATA
      "ows_enable_request" "*"
      "wms_srs" "EPSG:4326 EPSG:3857 EPSG:2154 EPSG:310642901 EPSG:4171 EPSG:310024802 EPSG:310915814 EPSG:310486805 EPSG:310702807 EPSG:310700806 EPSG:310547809 EPSG:310706808 EPSG:310642810 EPSG:310642801 EPSG:310642812 EPSG:310032811 EPSG:310642813 EPSG:2986"
      "labelcache_map_edge_buffer" "-10"
      "wms_title" "osm france"
   END
   IMAGEPATH "/tmp/ms_tmp/"
   IMAGEURL "/ms_tmp/"
END
DEBUG 1
CONFIG "MS_ERRORFILE" "stderr"
CONFIG "PROJ_LIB" "/home/mwh/osm-demo/basemaps"
PROJECTION
   "init=epsg:3857"
END
LAYER
   TYPE POLYGON
   STATUS ON
   GROUP "default"
   DEBUG 1
   NAME "land0"
   PROCESSING "APPROXIMATION_SCALE=full"
   PROJECTION
       "init=epsg:4326"
   END
   MINSCALEDENOM 332808204
   MAXSCALEDENOM 99999999999
   DATA "data/TM_WORLD_BORDERS-0.3.shp"
   CLASS
       STYLE
           COLOR 242 239 233
       END
   END
END
LAYER
   GROUP "default"
   STATUS ON
   MAXSCALEDENOM 99999999999
   MINSCALEDENOM 332808204
   TYPE LINE
   NAME "borders0"
   PROJECTION
     "init=epsg:4326"
   END
   DATA "data/boundaries.shp"
   CLASS
      STYLE
         WIDTH 0.5
         COLOR "#CDCBC6"
      END
   END
END

I followed the Tutorial given in MapProxy website to start my  proxy by starting the proxy. The tutorial example works fine. But when I use my own .map
file it shows an error "could not get any sources: no image returned from source WMS"

My mapproxy.yaml is shown below

services:
  demo:
  wms:
    md:
      title: MapProxy WMS Proxy
      abstract: This is the fantastic MapProxy.
      online_resource: http://mapproxy.org/
      contact:
        person: Your Name Here
        position: Technical Director
        organization:
        address: Fakestreet 123
        city: Somewhere
        postcode: 12345
        country: Germany
        phone: +49(0)000-000000-0
        fax: +49(0)000-000000-0
        email: info@omniscale.de
      access_constraints:
        This service is intended for private and
        evaluation use only. The data is licensed
        as Creative Commons Attribution-Share Alike 2.0
        (http://creativecommons.org/licenses/by-sa/2.0/)
      fees: 'None'

sources:
  test_wms:
    type: wms
    req:
      url: http://osm.omniscale.net/proxy/service?
      layers: osm

  mapserver:
    type: wms
    supported_srs: ['epsg:3857', 'epsg:4326']
    req:
      url: http://localhost/cgi-bin/mapserv?
      layers: default
      map: /home/mwh/osm-demo/basemap/osm-google.map
      transparent: true
    coverage:
      bbox: -20000000.0,-20000000.0,20000000.0,20000000.0
      bbox_srs: 'epsg:3857'

layers:
  - name: mapserver
    title: Mapserver Test Layer
    sources: [mapserver]

Currently I am stuck in this problem please help me to find some answer.

Thank you.


texts:
I have a GeoExplorer 6000 series handheld unit. I have a Terrasync 2.52 Serial key and installation disk.

I installed Terrasync 2.52 on the handheld unit but upon trying to connect to the receiver for 20 seconds, it says "GPS Disconnected" and never shows it connects to a satellite (COM 3 integrated GNSS Receiver is selected, all other settings default)

(also, under Setup, it says "GPS" instead of "GNSS" for some reason)

Arcpad uses the integrated receiver on COM 3 perfectly without issue.

I have also installed terrasync 2.53 with the exact issue. Any ideas?


texts: Older British Geological Survey geology maps of Kenya; need complete projection information

texts:
I have several geologic maps of Kenya produced by the BGS circa 1993. Noted on each map is the following:

Projection:-Transverse Mercator
Spheriod:-Clarke 1880 (modified)

That's great, but that doesn't tell me exactly which projection to use in ArcGIS. I've tried these two:

Adindan_UTM_Zone_37N and
Arc_1960_UTM_Zone_37N

Both put the maps in sort of the right place. There's a shift of about 50m to the west and 300m to the north. More or less.

Does anyone have a clue as to what the full projection parameters might be for these maps?


texts: Export GeoTIFF with pixel min max value from GEE

texts:
this is a question folllowing up with my previous question Export raster with its min max range from GEE and Export raster with its min max range from GEE 2. I need to export the GeoTIFF with pixels in the raster layer keep their index value so I can use it in ArcGIS. Below is the code I tried from @Kuik's answer and his link:

var visualizeMinMaxNDVI = function(maskImage) {
  var minMax = maskImage.reduceRegion({
  reducer: ee.Reducer.minMax(),
  bestEffort: true,
});
  var visParams = {
          min: minMax.getNumber('NDVI_min'),
          max: minMax.getNumber('NDVI_max'),
          palette:'green'
          };
  return maskImage.visualize(visParams)
            .set({min: minMax.getNumber('NDVI_min'), max: minMax.getNumber('NDVI_max')});
};

// set visualization parameters based on min and max value
var NDVIvis = maskImage.map(function(maskImage){
  return visualizeMinMaxNDVI(maskImage);
});
print(NDVIvis);
// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: maskImage,
  description: '21072016UNDVI',
  scale: 10,
  region: table,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});

I got error maskImage.map is not a function.

Also, how can I put visParams into Export.image.toDrive now (since visParams is a parameter in the function).
If you want to see my previous code please check the links at the start. Sorry for codes all over the place....


texts: One-to-many feature class to table in ArcGIS MakeQueryTable

texts:
I try to put up a one-to-many relationship between a feature class (one) and a table (many) using MakeQueryTable tool. One entry in the feature class needs to be linked to many entries in the table and I tried to set this up using SQL in the tool (I followed this ESRI article: http://support.esri.com/en/knowledgebase/techarticles/detail/37544):

[Feature_Class].FIELDA = [Table].FIELDB`

The records are matching, it just needs to relate them.

However, it gives me the following error:

  could not be retrieved properly. [Table.FIELDB] Failed to execute
  (MakeQueryTable).

If I try it the other way around, i.e.

[Feature_Class].FIELDB = [Table].FIELDA

it gives

  could not be retrieved properly. [Table.FIELDA] Failed to execute
  (MakeQueryTable).

As far as I can tell, everything is in perfect order, i.e. the table and feature class are in the same DB and I added the feature class first to the tool, in the hope to get out a feature class.

Anyone has any idea?


texts:
I'm working with ArcGIS and the following image is showing the polygons that I have in my shapefile. The seleted polygon should actually be this selected one + the blank area right below it.
This blank area has no attribute related to it. I mean, it isn't a polygon. It is just a hole in my shapefile.
Eventually, how do I use this shape to be part of my selected polygon?
I forgot to mention, but I have many holes like this in my shapefile. I'm wondering if there is a faster way of filling these holes than digitizing all those individual polygons!


texts: PostGIS or QGIS: Convert unsorted points to a single line by connecting each 2 closest points

texts:
I have shapefile/PostGIS table of points.

I want to connect them to form a single line.

When trying with PostGIS:

Select ST_MakeLine(points_table.geom) FROM points_table;

The problem is that because the points are not ordered, the line is crossing itself.

When trying in QGIS Points2One plugin - same issue.

After looking at the data - It seems that if connect each 2 closest points - the desired line will created.


texts: Why are labels and looks different between print layout and map canvas at same scale?

texts:
I know that some of my items uses map unit as size unit, others mm and labels points. But the map canvas where I do all my changes should still look exactly the same in print layout if I'm using the same map scale. Am I overlooking something?

I'm running QGIS 3.4.4 on macOS 10.14.3


texts: Have DEM, need polygon outlining region above a certain elevation

texts:
I have a DEM in HGT format. I need to isolate all areas above a certain elevation and view them as polygons. Something like this (let's say the pink polygon is everything above 1000m):

Seems like this would be a common GIS task, but all I found is the GRASS r.lake function. What I need is the opposite of a lake. :) QGIS/Python/GRASS solutions are preferred.


texts: Calculating area of overlapping polygons using ArcGIS for Desktop?

texts: In ArcGIS, I have a shapefile for a pipeline and another for plots of land that the pipeline shapefile touches and I was wondering if there is an obvious way to calculate the area that the pipeline covers within each of the individual polygons for the land plots shapefile.


texts: MODIS burn area - what do the pixels actually show?

texts:
I am using the MCD64A1.006 MODIS Burned Area Monthly Global 500m dataset to look at wildfire history over a certain area. I have downloaded the images from GEE and then uploaded them into Arc desktop. Within the images which display fire events, the pixels classified as burnt vary in their values. I am confused about what these values show. In some things I have read, it seems as though these represent the day of the fire, but the images are all for the first of each month so I don't understand how another temporal aspect really fits in.
If anybody has used this imagery and could explain to me what the varying pixel values mean that would be extremely helpful!


texts: Python comtypes and ArcObjects: Snippets GetApp() fails at 10.2

texts:
I have been experimenting with Python's comtypes module to access ArcObjects.  We are still using ArcGIS 10.0 as we are in the process of upgrading our ArcGIS for Server to version 10.2.2.  I am currently running ArcGIS desktop version 10.0 and have a VM with 10.2 installed for testing.

Anyways, the following from the Snippets module works to access an open mxd in version 10.0, but fails in 10.2:

def GetApp(app="ArcMap"):
    import arcpy
    """app must be 'ArcMap' (default) or 'ArcCatalog'\n\
    Execute GetDesktopModules() first"""
    if not (app == "ArcMap" or app == "ArcCatalog"):
        print "app must be 'ArcMap' or 'ArcCatalog'"
        return None
    import comtypes.gen.esriFramework as esriFramework
    import comtypes.gen.esriArcMapUI as esriArcMapUI
    import comtypes.gen.esriCatalogUI as esriCatalogUI
    pAppROT = NewObj(esriFramework.AppROT, esriFramework.IAppROT)
    iCount = pAppROT.Count
    if iCount == 0:
        return None
    for i in range(iCount):
        pApp = pAppROT.Item(i)
        print pApp
        if app == "ArcCatalog":
            if CType(pApp, esriCatalogUI.IGxApplication):
                return pApp
            continue
        if CType(pApp, esriArcMapUI.IMxApplication):
            return pApp
    return None

In 10.2, I get the following error message:

Traceback (most recent call last):
  File "\\arcserver1\GIS\_Resources\ESRI\Python\BMI_Library\bmi_arcobjects\_snippets.py", line 264, in <module>
    GetApp()
  File "\\arcserver1\GIS\_Resources\ESRI\Python\BMI_Library\bmi_arcobjects\_snippets.py", line 86, in GetApp
    pAppROT = NewObj(esriFramework.AppROT, esriFramework.IAppROT)
  File "\\arcserver1\GIS\_Resources\ESRI\Python\BMI_Library\bmi_arcobjects\_snippets.py", line 128, in NewObj
    ptr = CreateObject(COMClass, interface=COMInterface)
  File "\\arcserver1\GIS\_Resources\ESRI\Python\BMI_Library\comtypes\client\__init__.py", line 235, in CreateObject
    obj = comtypes.CoCreateInstance(clsid, clsctx=clsctx, interface=interface)
  File "\\arcserver1\GIS\_Resources\ESRI\Python\BMI_Library\comtypes\__init__.py", line 1145, in CoCreateInstance
    _ole32.CoCreateInstance(byref(clsid), punkouter, clsctx, byref(iid), byref(p))
  File "_ctypes/callproc.c", line 936, in GetResult
WindowsError: [Error -2147221164] Class not registered

I thought this may be an issue with versions since in the gen folder all wrappers ended with "_10_0".  Everyone's sys.path has access to the same comtypes which is on our server, so I blew away all the modules in the gen folder and called the following function to load all esri modules to ensure I was getting the most up to date wrappers:

def load_all():
    '''loads all object libraries'''
    from comtypes.client import GetModule
    mods = glob.glob(os.path.join(GetLibPath(), '*.olb'))
    for mod in mods:
        GetModule(mod)
    return

This loaded everything, and I still get that same error when running at 10.2.  I saw this post with a similar problem but the solution did not work for me: Python, comtypes and ArcObjects: Error creating AppROT object

Has anyone else had this problem?  If so, how did you alter the GetApp() function. Thanks.


texts:
I'm curious if it is possible to toggle on and off 'draft mode' in layout view in ArcGIS using python scripting.

Does anyone have any ideas?

Python is preferred, but if there is another solution in a different language it would also be beneficial!


texts: QGIS canvas refresh & rendering control from plugin?

texts:
Plugin Description

I have been working on a plugin which essentially doing the same  as "GPS information" live tracking but with capability of displaying two GPS streams. I am comparing streams and doing some other calculations on the run.

Feature

I am trying to find a way to play with rendering options. If you guys can recall there is a render check box at the bottom of QGIS, which can be used to disable rendering, I want a similar control from plugin but in my case I want to render the complete Layer Extent (not visible extent) before disabling the rendering option.

GOAL

I have set "Scale dependent visibility" on my layers and everytime I zoom in or out, I need to refresh the canvas to cater for new Layers which might appear or disappear at different zoom level. The issue is I want to disable rendering inside a particular zoom range after having rendered the whole layer extent. So what I need is a way to render to complete layer extent? How can I do that?.  If I don't have rendering to layer extent, my canvas looks like this when GPS moves out of rendered image:

To stop/start rendering I am using:

self.canvas.freeze(True)

I have tried following to get full extent of the layer:

ms = self.canvas.mapSettings()
extent_of_layers = ms.fullExtent()
ms.setExtent(extent_of_layers)

Haven't made any progress after above. any pointers will be helpfull.

Reply to First comment
I may not be accustomed to terms used here..sorry for that.. In short I want to know if its possible to take a snap shot/render a complete layers extent so that when my path/poniters move on the canvas, it will always have some background available for it, to have an idea where GPSs are on the map. Also to add, I am currently handling it by calling canvas.refresh() whenever pointer/path moves outside the visible canvas But it cause white flashes which is not easy on the eye especially zoomed in. So I was  thinking of rendering everything on Layers (even outside the visible canvas) at particular zoom levels and disable rendering otherwise, so that even when marker moves away from current extent there always some background available.


texts: stack raster using GDAL and numpy dstack

texts:
I am testing Savitzky-Golay filter with my MODIS dataset. Currently I am testing for 10 rasters (but later on will use my whole dataset) only just to test how should I create my script. So my initial workflow is to:

Access my raster files
each file I need to convert to 2D array using GDAL or GDAL Numeric
Append all those created 2D arrays to a list
Implement numpy.dstack() to stack my 2D arrays and hopefully I can begin some analysis.

However, I getting error in the 4th step (it says memory error) -- converting my list of 2D arrays to a stack. Could you guide me some advice how to proceed. I am just building my script so for now its not so much.

import os, sys
import gdal, numpy as np
import gdalnumeric as gd
from scipy.signal import savgol_filter

rasters = list()

ws = 'G:/Test/Raster'
for folder, subs, files in os.walk(ws):
    for filename in files:
        aSrc = gd.LoadFile(os.path.join(folder,filename))
        rasters.append(aSrc)

stackRast = np.dstack(rasters)


texts: Connecting two points with same label using QGIS?

texts:
I'm doing an MRR (mark-release-recapture) study on an insect. When I collect the insect I take the GPS point and I mark it with a unique number. When I recapture it I'm able to recognize it and I take the second GPS point to know how far it is from the first point. The column has the code of the "name" of the insect collected, for example N1.

I would like to connect the point with the same code N1 (The point where I collect the individual for the first time and all subsequent points where it was captured) to calculate the total distance travelled.