texts: Criterion for measuring the distance between a group of points to origin?

texts:
I wonder if is there any criteria to represent the distance of a group of points (in 2D space) to origin?

A simple criteria is the mean of distance of each point. But as can be seen in the picture, this criteria cannot capture the pattern of spatial distribution of points. Another criteria may be the square root of sum of squares of distance of points. But this is my idea and I cannot find any acceptable criteria in academic documents.


texts: Safety of using GDAL OpenFileGDB (reverse engineered) driver?

texts:
http://www.gdal.org/drv_openfilegdb.html

My understanding is that GDAL offers an alternative to accessing Geodatabase files via Esri's API SDK. Apparently Even Rouault reverse engineered the Geodatabase format, at least enough to get the basic geometry out (https://github.com/rouault/dump_gdbtable/wiki/FGDB-Spec).

It seems risky to use this (a) because perhaps the reverse engineering is imperfect (there are many "unknown"s in that spec) and maybe it could mess up on some unusual file because Even didn't understand the purpose of some field/flag, and (b) because there could be changes to the Geodatabase format at some point that the official SDK deals with correctly but the reverse engineering messes up on.

However, both of those seem fairly unlikely. And the big advantage for me of avoiding the Esri API SDK is that I don't have to integrate the libraries into my code, which is a big deal on Mac because the SDK only appears to be available in 64-bit, so will force my whole Mac build to become 64-bit, which is not something I want to deal with right now.

Does anybody have experience using the reverse engineered solution?


texts: How can I convert data in the form of lat, lon, value into a raster file using R?

texts:
I have a data set of values over a km grid in the continental U.S. The columns are "latitude", "longitude", and "observation", e.g.:

"lat"    "lon"     "yield"
 25.567  -120.347  3.6
 25.832  -120.400  2.6
 26.097  -120.454  3.4
 26.363  -120.508  3.1
 26.630  -120.562  4.4

or, as an R data frame:

mydata <- structure(list(lat = c(25.567, 25.832, 26.097, 26.363, 26.63),
lon = c(-120.347, -120.4, -120.454, -120.508, -120.562),
yield = c(3.6, 2.6, 3.4, 3.1, 4.4)), .Names = c("lat",
"lon", "yield"), class = "data.frame", row.names = c(NA, -5L))

(the full data set can be downloaded as csv here)

The data are output from a crop model (intended to be on) a 30km x 30km grid (from Miguez et al 2012).

How can I convert these to a raster file with GIS - related metadata such as map projection?

Ideally the file would be a text (ASCII?) file because I would like for it to be platform and software independent.


texts: Using Tilestache for serving maps to iPAD

texts:
I am planning to use tilestache as the tile server (map server) for serving maps to both iPAD and web clients.

I plan to generate tiles by Tilemills and store tiles in a Sqlite database (mbtiles) and serve it using Tilestache.

My questions are:

(A) How scalable is the solution, for example the project may have 200 IPAD users in field?

(B) Is client side local caching supported or do we have to pre-load tiles in the local storage?

(C) Most important, can any such mobility solution using tilestache be sited to see the performance results?


texts: How to convert ASCII file to TIN where XY is YX in ArcGIS?

texts:
I have some ASCII files (you can download it from: http://analizasan.kvprojekty.pl/underwork/gisse/M3451Ca1.zip) that I would like to create TIN.

First I used "ASCII3DToFeatureClass" tool to create shapefiles with points. Then "CreateTin" tool to create TIN. Everything is fine but my XY coordinates are reversed and my TIN is somewhere else compared to the real.

Can I indicate X and Y when I create TIN by "CreateTin" tool? I know that I can edit .asc file in TXT Editor but it is a lot of work.

These files are in PUWG_1992 coordinate system (Poland):
Projection: Gauss_Kruger
False_Easting: 500000,0
False_Northing: -5300000,0
Central_Meridian: 19,0
Scale_Factor: 0,9993
Latitude_Of_Origin: 0,0
Linear Unit: Meter (1,0)

Geographic Coordinate System: GCS_ETRS_1989
Angular Unit: Degree (0,0174532925199433)
Prime Meridian: Greenwich (0,0)
Datum: D_ETRS_1989
  Spheroid: GRS_1980
    Semimajor Axis: 6378137,0
    Semiminor Axis: 6356752,314140356
    Inverse Flattening: 298,257222101

I'm waiting for help.


texts: Dividing range of territory into squares with different colors from points with accident information using QGIS

texts:
I have points with accident information and I need to divide the range of my territory into squares with different colors.
And the coloring of the square must depend on how many points (accidents) lie within its boundaries.
How can I achieve this using QGIS?
The picture above is an example of what I want to make


texts: Writing numpy array to rasterio raster file generates empty file

texts:
I'm using Rasterio to read GeoTIFF from Sentinel 2 and calculating NDVI, but when I try to save as a new TIF file, the file turns to transparent and empty file. What am I missing?

That's the code:

class MapService:

    def __init__(self):

        dsRed = rasterio.open('data/bands/B04.TIF')
        bandRed = dsRed.read(1)

        dsNIR = rasterio.open('data/bands/B08.TIF')
        bandNIR = dsNIR.read(1)

        ndvi = np.zeros(dsRed.shape, dtype=rasterio.float32)

        ndvi = (bandNIR.astype(float)-bandRed.astype(float))/(bandNIR+bandRed)

        kwargs = dsRed.meta
        kwargs.update(
            dtype=rasterio.float32,
            count=1,
            compress='lzw')

        self.ndvi = ndvi
        with rasterio.open('result/example-total.tif', 'w', **kwargs) as dst:
            dst.write(ndvi.astype(rasterio.float32),1)


texts: Automatic text wrapping in ArcMap layout text element?

texts:
How can I word-wrap a text element on an ArcMap layout, when it contains a lot of text?

go into Layout mode
choose > Insert > Text
paste in a long text block eg lorem ipsum

The text block now stretches across the screen, and there's no apparent way to change the width of the text block. Is there a way to set the text block to the desired width, and have ArcMap automatically enter line breaks?

(The way that, um, every other piece of software which handles text does...)

The workaround is to enter line breaks manually, which is great until you decide to change the width, font size, etc.


texts:
I'm using osm2po to break up OSM map into routable segments, and pretty much always I want the streets to go from one intersection to another (or till some alley).

However, if underlying OSM map has a street that is already broken up into multiple segments (a lot of times arbitrarily), I'm ending up with the same segments in my map, so even some streets with no alleys coming into them are split up in two or three parts. Is there some way to instruct osm2po to merge these smaller segments into one?

Example:

Map that I generate, with Kolpingplatz separated into two segments (orange and red).
http://imgur.com/a/w1p0o

Same street in OSM, composed of two segments:

http://imgur.com/a/txsVH

And in reality, I just want that segment of Kolpingplatz to be one street segment, not two.


texts: Calculating Zonal Statistics outside of polygon using ArcGIS Desktop?

texts:
Is there a way to calculate zonal statistics outside of polygon?

I want to extract the mean values outside of the polygon. In this case, instead of calculating for the land, I want to use the ocean.

I use ArcMap 10.5 for my analysis.


texts: QGIS Python, adding attributes and populating the fields

texts:
I am writing a script to add attributes and then populate the fields. When I run the script in QGIS in the python console, it works. But when I put it in a script and run it through a plugin, it only adds the fields, but doesn't populate them.

Here is the script,

def addAtts(attName, varType):
    layer = iface.activeLayer()
    layer.dataProvider().addAttributes([QgsField(attName, varType)])
    layer.updateFields()

attDict = {'COLOR': 176, 'ICON': None, 'L_LIMIT': 0, 'LAYER': 50, 'MAJ_CAT': 'Nautical', 'MIN_CAT': '(AREA)Lake', 'MODE':  'N', 'VALUE': None, 'U_LIMIT': 34, 'FONT': '14B', 'P_LIMIT': 0, 'P_ICON': None, 'BFR_COLOR':  0, 'BFR_WIDTH':  0, 'RTP_LIMIT': 34, 'RTP_COLOR': 43, 'WIDTH': 1, 'BDR_COLOR': 43, 'BDR_WIDTH': 1, 'IMG_RENDER': 255, 'BFR_PATTERN': None, 'PATTERN': None, 'BDR_PATTERN': None, 'FILT_SRCH': None }

layer = iface.activeLayer()
layer.startEditing()

num = 0

for x in attDict.keys():
    num += 1
    it_is = str(attDict[x])
    try:
        int(it_is)
        it_is = True
        #qType = 'QVariant.Int'
        addAtts(x, QVariant.Int)
        layer.changeAttributeValue(0,num,attDict[x])

    except ValueError:
        it_is = False
        #qType = 'QVariant.String'
        addAtts(x,QVariant.String)
        layer.changeAttributeValue(0,num,attDict[x])


texts:
Does anybody know if Scimap ( http://www.scimap.org.uk ) works with QGIS?

I know it is supposed to operate with ArcGIS and SAGA but as I am very much a beginner with QGIS I would like to know if its possible to run the software before I spend a lot of time trying to get it going.

I am running QGIS 2.0.1-Dufour on a rather elderly 8 Mb Macbook Pro


texts: How to set monitoring properties in order to moniitor a rest interface with PRTG?

texts:
I am trying to use PRTG as a monitoring tool for ArcGIS Server 10.3.1 Rest interfaces.

I am having problems finding the correct sensor to add to the PRTG Map in order to monitor a Rest endpoint correctly.

In the properties section on the sensor I have to enter 2 different values

URL and

XML - Node (and optional Properties).

I have entered the URL of a service as displayed in the picture.

I am a bit confused about the XML-Node - as the rest interface returns JSON.


texts: How to check-in ( release, make it free) ArcGIS10.2 desktop license using arcpy?

texts: There is a way to do this using ArcGIS license Server Administrator, but I need a command to use it in arcpy scripting.


texts: QGIS - “Error Entry Point Not Found” when launching QGIS 1.8

texts:
I am trying to download QGIS 1.8 on my windows operating system. When I tried to launch QGIS 1.8 a dialog box appears, preventing me from launching the software. This is the error that is generated:

I don’t know what this error message mean or how to fix it.


texts: Why do I get value of GDAL GetMaximum() different than that of numpy array max() for same image band?

texts:
Here is a code snippet

from osgeo import gdal
import numpy

file = gdal.Open('DEMO.TIF')
band = file.GetRasterBand(1)
max = band.GetMaximum()
if not min or not max:
    (min,max) = band.ComputeRasterMinMax(True)
print(max)

redArray = band.ReadAsArray()
maxArr = redArray.max()
print(maxArr)

Apparently both the max values need to be same as I am working on the same dataset and band. But I get the result as,

49820.0
53041

Why are these values not same?


texts: Combining values in shapefile using ArcMap?

texts:
I have a shapefile with 10 different values (geological ages) which I want to combine into four values (into broader age categories).

Is there a tool which I can use to selectively combine values in ArcMap?


texts:
I need to filter a Sentinel-2 image collection in order to have almost no clouds.
I used the following code:
//Dates of interest
var START = ee.Date('2017-03-28');
var FINISH = ee.Date('2020-07-31');
//Merge different plot
var Plot = Plot01.merge(Plot02);

//Create ImageCollection
var ImageCollection = ee.ImageCollection(S2L2A)
.filterDate(START,FINISH)
//Filter data by location
.filterBounds(Plot)
//Filter data based on metadata property -> cloudy pixel percentage <5%
.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 5)
//Return clipped image
.map(function(image){return image.clip(Plot)})

The filter of cloudy pixel percentage <5% refers to the entire Sentinel-2 tile. How can I perform this filtering process referring just to my Plot?


texts: "Definition Query on a versioned Geodatabase"

texts:
I would like to implement Specific views to the SDE users,  based on Definition Queries of common attributes , like Date created > = October 1.

But on first look, the definition queries are for "feature Display setting" of each featureclass - Our aim is to implement it for the Feature dataset in a particular version

The below article did not help, as I could not locate , choose Production Cartography In the Categories list.
http://resources.arcgis.com/en/help/main/10.1/index.html#//0103000000s5000000

What is the procedure to display features on a common attribute, across the feature dataset in a Version.


texts: How it can be made a spherical cap surface with qgis raster calculator

texts:
I would like to make a test in QGIS raster calculation: I'd like to make in a flat raster surface (0m as elevation) a kind of hill or dome with the spherical cap geometry formula.

P is the point vertex projection of the height in the base of the cap, the coordinates of which are x, y.

Could anyone help me?


texts: AttributeError: OracleDBConnector instance has no attribute 'cache_connection' in QGIS DBManager

texts:
I have a connection with an Oracle Database but when I try to make a query  using DB Manager of QGIS it appears "Not Connected" and I get a Python error:

1   Traceback (most recent call last):
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_model.py", line 426, in rowCount
                self._refreshIndex(parent, True)
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_model.py", line 477, in _refreshIndex
                if item.populate():
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_model.py", line 158, in populate
                if not connection.connect():
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_plugins\oracle\plugin.py", line 126, in connect
                return self.connectToUri(uri)
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_plugins\oracle\plugin.py", line 73, in connectToUri
                self.db = self.databasesFactory(self, uri)
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_plugins\oracle\plugin.py", line 79, in databasesFactory
                return ORDatabase(connection, uri)
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_plugins\oracle\plugin.py", line 161, in __init__
                Database.__init__(self, connection, uri)
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_plugins\plugin.py", line 211, in __init__
                self.connector = self.connectorsFactory(uri)
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_plugins\oracle\plugin.py", line 164, in connectorsFactory
                return OracleDBConnector(uri, self.connName)
              File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\db_manager\db_plugins\oracle\connector.py", line 112, in __init__
                if self.cache_connection:
            AttributeError: OracleDBConnector instance has no attribute 'cache_connection'

It is strange because I can add the layers using QGIS Browser, so I guess it is connected.

Anyone knows how to solve it?


texts: Mouse single click problems after QGIS install

texts:
Has anyone had an issue with their single click becoming a double click after installing QGIS on Dell machines?

On two Dell workstations I have installed QGIS and on both my mouse has a tendency to act as a double click unless I hold down on the left button much longer than normally required.

I have adjusted the mouse click settings with no affect.

I have narrowed it down to QGIS and it does not get resolved until I do a total clean install of Windows. Removing QGIS from the machine does not remedy the issue--only a clean install of the OS works. The issue affects all mouse clicking in all programs to include file explorer.

I run a Windows 10 Home 64 bit Dell Inspiron 3847 i7-4790 @ 3.60ghz w/16GB ram. I have tried switching to a mouse with its own driver trying to over ride QGIS and doing a pure windows 10 install with no Dell components.

This issue does not occur on my Toshiba laptop.


texts: PostGIS: Interpolate/Cluster geographic points to average points

texts:
Setup: PostGIS 3.1, QGIS 3.16
I have a bunch of geographical data (geometries). They show more or less the same route, just measured several times. Due to the many measurements and the GPS precision, the measured points differ, of course:
General visualization of the data:

A more detailed view:

Next to the coordinates the data points keep some additional data, for example the speed of the vehicle which delivered the measurements. This is displayed by the color.
I am interested in the average speed at the track. So I did some interpolation (clustering?) and aggregated the points using PostGIS' ST_SnapToGrid() function:
SELECT
    ST_SnapToGrid(geom, 0.001) as c_geom,
    avg(speed) as avg_speed
FROM my_data
GROUP BY 1

This, naturally, yield a grid-like aggregation:

However, I am searching for a function which I can use like the ST_SnapToGrid() function, but which do not simply trunc the digits of the latitude and longitude but does more like an averaging, so the resulting points are more located on the track itself ("snap to track"). Of course, the speed value needs to be averaged as well.
Desired output:
I am searching for a function which averages lat/lon/speed values within a 10 meters radius or something like that:

I am not even sure if this is possible.


texts: Feature Layer not showing up when adding an automatic rotation function to SceneView

texts:
I am building a simple web application using the ArcGIS API for JavaScript (4.x). I created a map with an associated FeatureLayer and then a SceneView and everything works fine. However, I'm trying to make the scene rotate automatically. When I implement a function using scheduling the rotation works but the FeatureLayer dissapears and only the basemap is visible. Has anyone experienced the same problem or does someone know an approach to this issue?

Here is my code. The function that makes rotate the scene is at the end, inside scheduling.addFrameTask. I took it from an example in http://proceedings.esri.com/library/userconf/devsummit16/papers/dev_int_204.pdf (page 21).

<!DOCTYPE html>

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

    <link rel="stylesheet" href="https://js.arcgis.com/4.8/esri/css/main.css">


texts: How can I change layers visibility while I zoom in GeoServer?

texts:
I want to change layers visibility while I zoom , here are my screenshots.

where is my error ?


texts: How to return individual address components (city, state, etc.) from GeoPy geocoder?

texts:
I'm using GeoPy to geocode addresses to lat,lng. I would also like to extract the itemized address components (street, city, state, zip) for each address.

GeoPy returns a string with the address -- but I can't find a reliable way to separate each component. For example:

{street: '123 Main Street', city: 'Los Angeles', state: 'CA', zip: 90034, country: 'USA'}

The Google geocoding API does return these individual components...
is there a way to get these from GeoPy? (or a different geocoding tool?)


texts: Pixels to lat,lon in Python from .shp, and .prj

texts:
I have the task to analyse ice types. I'm using Python to assist, no other programs is possible, unless I can export pure text dumps from said program, it's free and it can be fully automated and interfaced with Python.

Specifically, I'm using data from http://nsidc.org/data/g02171 (Canadian ice service), Western Arctic area. I understand the shapefile. The .prj file I don't fully understand. I realize it is parameter-information to show how the flattening / projection onto earth transformation is done. I've found http://docs.opengeospatial.org/is/12-063r5/12-063r5.html#49
but honestly, that is overwhelming.

What do I need to do to apply the inverse transformation?

I'm guessing I need information about the transformation, if so, where?


texts:
I am working with vector geometries that I recently simplified using the native:simplifygeometries method in the processing toolbox. I am pleased with the simplification but I want to "fill" the gaps in between the geometries so there is no empty space.
How do I make all of the geometries kiss each other?
I don't care which geometry the extra space is allocated to but it cannot itself be a new geometry.
Perhaps there is a better simplification method to use that doesn't lead to this problem, but I need the features to be highly simplified.


texts: Importing geocentric cartesian coordinates csv file in QGIS

texts:
I'd like to import a CSV files with geocentric (ITRF) X,Y,Z coordinates points, like

 CHPH , 4236233.08156, 110998.26463599999, 4751117.477176
 MAN2 , 4274275.799144, 11584.521544, 4718386.149148
 MLVL , 4201576.823332, 189860.284372, 4779064.903872
 SIRT , 4213550.772464001, 162494.69744, 4769661.8907079995
 SMNE , 4201791.9088, 177945.66576799998, 4779287.023676001

But the Text Layer Import assistant only allows 2D/projected coordinates

Any idea to how importing easily those points without preliminary conversion ?
(The destination is a OSM layer in EPSG3857 Web Mercator projection system)


texts: Why North direction is shown in red color on aspect maps?

texts:
I have been recently exposed to an aspect map that looked like this:

I wonder why North-y orientations (N, NE, NW) are shown in red and orange while South-y orientations (S, SE, SW) are shown in light and dark blue?

According to the color theory, blue is a cold color and red is a warm color. If I understand things right, in the Northern hemisphere, areas with the North aspect are exposed to the sun less that the ones with the South aspect. This makes me think of areas with South aspect as of "warmer" areas and ones with North orientation as of "colder" areas, so I would reverse the color schema used.

Is there any particular reason why these colors are used as a standard for aspect maps? I have browsed through many aspect maps online (link1, link2, link3) and all of them use the same color ramp I referred above.


texts: What Coordinate Reference System does iNaturalist export observations in?

texts:
iNaturalist has a nice export tool which selects up to 200,000 observations based on a query and writes them as a clean CSV, including latitude and longitudes where provided. To import that CSV file into QGIS, a Coordinate Reference System needs to be defined.
I cannot find what CRS iNaturalist uses to store lat/long. Since this CSV seems to import correctly into Google Earth with no problems, I presume iNaturalist uses the same CRS as Google (EPSG:900913). It would be good to confirm this but I haven't found the info yet so I am asking here.


texts: Difference between Save As and Convert in Surfer?

texts: Is there any difference in the file saved when using "Save As" compared to "Convert" to save/convert a Surfer .grd as an ASCII .asc grid?


texts: I updated a plugin such that now it does more operations than the original title suggests.  Is there an easy way to change the name without going through all of the files line by line?


texts: Calculating areas in metric units on polygons that are stored in geographic coordinate reference system using ArcPy?

texts:
If one uses arcpy to calculate areas for a field in metric units on a polygon shapefile in a geographic coordinate system it will work, and give you plausible results. According to the documentation it uses 'a geodesic algorithm' to calculate distances but the documentation doesn't talk about how areas are calculated. How is it producing results in square kilometers in the following code and how does this calculated area compare in accuracy compared to projecting the data into an equal-area projection before calculating the area?

result_file = settings.get_boundary_file_results_location(boundary_short_name)
area_field_name = settings.NEW_FIELD_NAMES['Area square kilometers']
arcpy.AddField_management(result_file,
                          area_field_name,
                          "FLOAT")
expression = "{0}".format("!SHAPE.area@SQUAREKILOMETERS!")
arcpy.CalculateField_management(result_file, area_field_name, expression, "PYTHON_9.3")

And the properties of the shapefile:

Edit: There is a closed but nonetheless useful question here that deals with some of this: ArcGIS Length and Area Calculation Scenarios


texts:
I'm a beginner at QGIS and I'm having trouble aligning my shapefiles.

I have two datasets: one originally in WGS84 (lat/long) and another in Corrego Alegre/UTM zone 24S.
I need to convert the second into WGS84 for my work, however if I go to Vector -> Data Management Tools -> Define Current Projection and choose to do so, I'm told this is an invalid coordinate system.

By going to Save as... I'm able to save a shapefile that seems to be in WGS84, but it is dislocated from the other polygons of my first dataset.

How can I get these files to align??
This is a screenshot of the mis-alignment, even though the novo-recifes-v2.shp (coral reefs) is set as WGS84, just like the continent (global_polygon.shp)


texts: Geowebcache automatic seed is not working

texts:
I created new geowebcache task via form as following:

But my task list is empty:

I want to create tiles automatically. How can I do this?

UPDATE:

I noticed that geoserver seed is seeding geoserver default sample data. For example I seed topp:tasmania_state that source is file. But I am creating postgresql publish from database, this is not working. EPSG is 4326. But gwc demo page is creating tiles.


texts: Getting download statistics from ArcGIS Hub

texts:
I have configured a new GIS page within ArcGIS Hub.  It contains links to other GIS sites as well as to local road closure apps, river levels, and other local information.  It also contains links to download our GIS data broken down into categories.
Is there a way to find out how many times our GIS data has been downloaded per feature within a given time period (i.e. downloaded centerlines 47 times in 4th qtr 2020)?


texts:
I'm currently doing a project on finding blocked roads in a road network layer. I've figured that I can use the topology checker with the rule 'must not have dangles' and find them this way. But I can't figure a way to export this to a layer or any other file format I can use for further analysis. Any ideas?

Please see the attached image of my result then it might be more clear what I mean (the QGIS version is in Danish, please let me know if it's a problem)

The topology-checker I'm refering to is the one in Vector -> Topology-checker where you have the option to define rules.

I should add that my road network is combined to 1-line feature and not the seperate road sections. In this way I'm able to find the errors with the 'must not have dangles' rule.

I am aware of ArcGIS has this export topology option, but since I'm mainly using QGIS and Mapinfo I prefer to find a solution through these.


texts:
I have two polygons. One contains fields(X,Y,Z) and the other contains soil types (A,B,C,D). I want to know what  area of every field contains which type of soil. I tried the following:

library(rgdal)
library(rgeos)
Field<-readOGR("./","Field")
Soil<-readOGR("./","Soil")
Results<-gIntersects(Soil,Field,byid=TRUE)
rownames(Results)<-Field@data$FieldName
colnames(Results)<-Soil@data$SoilType

> Results
      A     B     C     D
Z  TRUE FALSE FALSE FALSE
Y FALSE  TRUE  TRUE FALSE
X  TRUE  TRUE  TRUE  TRUE

and achieved good results with it telling me which field contains which soil type. However, how do I get the area instead?


texts:
There was a lot of svg-files included when installing QGIS (C:\Program Files\QGIS 3.4\apps\qgis-ltr\svg), but the program can't really seem to find these.
I've made sure the path is correct under Settings/Option/Systems and i've also tried adding them to the active profile folder (C:\Users\UserName\AppData\Roaming\QGIS\QGIS3\profiles\default\svg), and made sure the path there is correct as well.

When trying to import them using the Style Manager, i seem to only be able to choose xml-files.

How can I access the files which are included in the installation?


texts: Reading non-standard JSON / TopoJSON files in GeoPandas

texts:
I am unable to read a TopoJSON / JSON file in GeoPandas. The issue is that it's a JSON file with keys that are "probably" not standard.
I tried to use the from_features() method but it gives the error: KeyError: 'geometry'
Getting the GeoDataFrame is fine, but it doesn't have the geometry column as one would like.
File Source
File's data structure -
state:  "GA"
stateName:  "Georgia"
voteData:   […]
topojson:
    type:   "Topology"
    arcs:   […]
    transform:  {…}
    objects:
        counties:
            type:   "GeometryCollection"
            geometries: […]
        stateOutline:   {…}
countySummaryStats: {…}

What I tried:
import geopandas as gpd
import pandas as pd
import requests
uri = "https://static01.nyt.com/newsgraphics/2020/10/29/election-early-battlegrounds/030056173a56b8c767e453b62d8cee320984b92d/GA.json"
r = requests.get(uri)
data = r.json()['topojson']['objects']['counties']['geometries']
gpd.GeoDataFrame().from_features(data)

When trying gpd.GeoDataFrame().from_features(r.json()['topojson']['objects']['counties']) I get the error -
TypeError: string indices must be integers

Ideally I'd like to have voteData's fields corresponding to geometries.
Any suggestion how to read this file?


texts: Relating shapefile to excel table in ArcMap?

texts:
I have successfully joined a point shapefile to .CSV file in ArcMap 10.0, and the type of the fields which are related in the two files are the same(text, length:255). But when I remove join and relate the two files, open the attribute table of the shapefile, and click the "Related Table" button, the option is greyed out...
How can I solve this problem?


texts: Removing the underlying polygons to only get the top polygon

texts:
I'm trying to do a suitability analysis in ArcGIS Pro but I've seem to run into a problem with overlaying polygons that shouldn't be overlaying. I put all the data into a weighted sum and then divided the weighted Sum into 11 different "classes" using raster calculator. I should then get the areas of all classes perfectly alligned next to each other but they aren't (as shown in the foto's).

I'm trying to only get the top layer if it overlaps with the layers beneath but I can't seem to figure out how. Does anybody know how to do this?


texts: Download metadata .xml file from geonode / geoserver

texts:
I am new to both geonode and geoserver. I have a problem with downloading xml file from layer that I uploaded.
For example, I have a.shp, a.shx, a.prj, a.dbf, and a.xml. I uploaded them to geonode. I check in geoserver data directory, the files are exist.
But, when I tried to download the layer in shapefile zipped format, it doesn't include the a.xml file.

I know that geoserver uses WFS. But I can't find option to download the .xml file.

How to do it? Can I access directly to the directory (which is not so good), or is there a way to do it?


texts:
I'm trying to create a function that will be essentially taking various parameters perform pgrouting with Dijkstra's algorithm, according to the user input as parameters. So the I'm making a minimal example which is the following:

 CREATE OR REPLACE FUNCTION _a_route(source integer, target integer)
 RETURNS SETOF record AS
 $$
 BEGIN
 SELECT * FROM pgr_dijkstra('SELECT id, source, target, time AS cost FROM roads_noded', $1, $2, false);
 END;
 $$
 LANGUAGE plpgsql;

This should be taking source and target as input and run the algorithm.

But I'm always getting the same error:

ERROR:  query has no destination for result data
HINT:  If you want to discard the results of a SELECT, use PERFORM       instead.
CONTEXT:  PL/pgSQL function _a_route(integer,integer) line 3 at SQL statement

So what am I doing wrong here? SETOF record is supposed to be the equivalent of returning a table

Postgresql Version:   9.6.1
PostGIS  : 2.3.1


texts: How can I show in a report or otherwise both the field name and field alias for a feature?

texts: I am trying to create a report that will list both the field name and the field alias for every table we have in our geodatabases.  I know there is no existing tool that can do this and I have found several python scripts that can do things such as batch edit the alias field but I am just looking to display both fields for each table or a null value in the case of no alias shown. I am working in Arc 10.3.1.


texts: using multiple mapbox wmts baselayers with dififerent labels in QGIS composer

texts:
I have created two basemap layers in mapbox.

https://api.mapbox.com/styles/v1/akrystos/cjr1epzf905a32sqhkwvjl5i7.html?fresh=true&title=true&access_token=pk.eyJ1IjoiYWtyeXN0b3MiLCJhIjoiY2o0YWI5Ynh0MTJleTMzbGNiMmNycXJtMyJ9.usNtifL5jfkUdC0NgNvChQ#5.0/0.400000/37.850000/0
https://api.mapbox.com/styles/v1/akrystos/cjr1f15ar91b72rpcgrgtzhb9.html?fresh=true&title=true&access_token=pk.eyJ1IjoiYWtyeXN0b3MiLCJhIjoiY2o0YWI5Ynh0MTJleTMzbGNiMmNycXJtMyJ9.usNtifL5jfkUdC0NgNvChQ#5.0/0.400000/37.850000/0

I am using basemap 1 in a main map in group 1 with my data layer.
I am using basemap 2 in an inset map in group 2.

The two maps appear in the composer correctly and are refreshed/locked. But, once I export the map to pdf, the labels in both maps revert to whatever is in the dataframe view of arcgis.

Any thoughts?

As a workaround, I exported the inset as an image and added it to the composer as a static image.

QGIS version 2.18.27
windows 10 pro

Amy

Please see the images of the data view, composer, and exported pdf below.

The composer and data view of the two maps

The composer and exported PDF of the two maps


texts: Dissolving while retaining field informations in the same entry instead of multiple entries

texts:
I have 1520 covenants. Many covenants have more than 1 VM symbols attached to it, currently showing as rows. EG. Covenant #97 has 4 VM Symbols, with the VM Symbols 1, 67, 117 and 171. How do I change the table to show Covenant #97 as a single entry with 4 columns, each column showing 1, 67, 117 and 171. There will be  entries with more or less columns populated. This is fine as those columns will be left as null. There will be a maximum of 8 VM Symbol Columns as 8 is the maximum frequency which occurs for a covenant. IS there a quick way of doing this? Otherwise, a model or a python script would also be welcome.

I thought if there is a way to manually add new fields and populate the field with the relevant symbols, I can do a dissolve at the end for FID Covenants so we don't end up with multiple similar entries.


texts: Where does the Geoid information in a GPS unit come from?

texts:
In the NMEA GGA-message the 11th field is Geoid separation. How does the GPS-unit  "know" this separation?

Are GPS units preloaded with an EGM (Earth Gravitational Model)?


texts: qgis: generalisation & covex - how to? what plugin?

texts: I'm looking for a way to generalize and simlify polygon to remove "internal spikes" as ilustrated below. The violet in input and green expected poligon.


texts: Using Data Driven Pages to clip to data driven page extent without clipping points at border of extent in ArcMap

texts:
I am creating a mapbook in ArcMap 10.6.1 using Data Driven Pages, and excluding data outside of the data driven feature specified. However, the points within the feature of interest always get clipped at the border of the polygon. I.e. I often get cut-off points as seen in the image below (the image represents a polygon with random points generated along its edge):

My ultimate question is how do I keep the point symbols from being cut off at the borders?

I should also note that I have three data frames: one inset, and two side-by-side maps. So the solution would need to affect at least two of the data frames.

This section is for those that need further explanation of how I got here:

I first turned on data driven pages using the Data Driven Pages toolbar and set up my data accordingly. I then went to View > Data Frame Properties > Data Frame tab > Clip Options box and changed the settings to Clip to current data driven page extent. I then clicked Exclude Layers... checking the layers I wanted to serve in my background. I did not check my points layer.

If you need any further details, feel free to ask.


texts: Add raster value to features AttributeError: 'module' object has no attribute 'TYPE_RASTER'?

texts:
When I try to run the "add raster value to features" (available into SAGA Tools), I get this error:

Traceback (most recent call last):
  File "/usr/share/qgis/python/plugins/processing/gui/AlgorithmDialog.py", line 249, in accept
    msg = self.alg._checkParameterValuesBeforeExecuting()
  File "/usr/share/qgis/python/plugins/processing/core/GeoAlgorithm.py", line 239, in _checkParameterValuesBeforeExecuting
    return self.checkParameterValuesBeforeExecuting()
  File "/usr/share/qgis/python/plugins/processing/algs/saga/SagaAlgorithm.py", line 343, in checkParameterValuesBeforeExecuting
    param.datatype == dataobjects.TYPE_RASTER):
AttributeError: 'module' object has no attribute 'TYPE_RASTER'

I'm using QGIS 2.18.10 and Python 2.7.12 on Ubuntu 16.04

I don't have any clue to solve this error.


texts: How can I get QGIS Plugins to work through the proxy internet settings?

texts:
For QGIS on my work computer, I add in the appropriate information under Settings - Network - Proxy.  This allows me to access plugin repositories, but sometimes plugins do not work.  I would get 'network failure' or proxy restrictions.

How do QGIS plugins that require internet access work if the Proxy/network settings fail?


texts: how to remove cloud values of Landsat image from R

texts: How can I remove cloud value of the Landsat8? How can I estimate cloud values  for eg: modis cloud values are from 700 to 800. By using R how can I estimate>


texts: GDAL for C# does not contain all methods as in C++

texts:
According to GDAL's official site, GDAL 2.1 OGRDataSource class is deprecated and GDALDataset should be used. Methods like OGR_DS_GetLayerByName should be replaced by GDALDatasetGetLayerByName().

Now in GDAL for C#, latest 2.1.0 version, I cannot find this new methods in the Gdal.Dataset class.

Where did they go?


texts: Finding lowest elevation point over borders of Polygon

texts:
I am new to QGIS and I have a polygon shapefile and a raster DEM file. I want to find the lowest elevation point over boundaries of each polygon in QGIS.

Can anyone please tell me the way to do this?

I also want to know the list polygons this points is part of (as it is boundary point). I want longitude an latitude of lowest elevation point over boundary of each polygon and as this point lies on shared boundary hence I need list of other polygons on which this point lies.

I need this to make a stream flow sequence for such drainage points and need to use this for other computation after having such dependencies of water stream flow among polygons

As in image encircled red point lowest elevation point for polygon numbed 25 and blue encircled point is for polygon no 18. Also Red encircled point lies in boundaries of polygon 5,7,25,18. I want similar lowest elevation points (Long. and Lat.) for each polygon and polygon list for each such point.


texts: Standard Deviational Ellipse with Open Source Python (GDAL/OGR etc)

texts:
I am trying to mimic the Standard Deviational Ellipse tool from ArcGIS through Open Source Python avenues. The formula is here.

At the moment I am failing to get the correct X and Y axis lengths for the ellipse, the first step in the formula.

When I run the tool in ArcGIS the values are 4607.796039, 7816.667667 but I get 3948.758414, 5057.00957766.

Here's my code so far...

from osgeo import ogr
from shapely.geometry import MultiLineString
from shapely import wkt
import numpy as np
import sys, math

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
## path to the FileGDB
gdb = r"C:\Users\******\Documents\ArcGIS\Default.gdb"
## ope the GDB in write mode (1)
ds = driver.Open(gdb, 1)

input_lyr_name = "Birmingham_Burglaries_2016"

output_fc = input_lyr_name + "_standard_ellipse"

## reference the layer using the layers name
if input_lyr_name in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
    lyr = ds.GetLayerByName(input_lyr_name)
    print "{0} found in {1}".format(input_lyr_name, gdb)

if output_fc in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
    ds.DeleteLayer(output_fc)
    print "Deleting: {0}".format(output_fc)

## for each point in the layer
## get the x and y value
## and place in an array
try:
    first_feat = lyr.GetFeature(1)
    if first_feat.geometry().GetGeometryName() in ["POINT", "MULTIPOINT", "POLYGON", "MULTIPOLYGON"]:
        xy_arr = np.ndarray((len(lyr), 2), dtype=np.float)
        for i, pt in enumerate(lyr):
            ft_geom = pt.geometry()
            xy_arr[i] = (ft_geom.Centroid().GetX(), ft_geom.Centroid().GetY())

    ## for lines we get the midpoint of a line
    elif first_feat.geometry().GetGeometryName() in ["LINESTRING", "MULTILINESTRING"]:
        xy_arr = np.ndarray((len(lyr), 2), dtype=np.float)
        for i, ln in enumerate(lyr):
            line_geom = ln.geometry().ExportToWkt()
            shapely_line = MultiLineString(wkt.loads(line_geom))
            midpoint = shapely_line.interpolate(shapely_line.length/2)
            xy_arr[i] = (midpoint.x, midpoint.y)

except Exception:
    print "Unknown geometry for {}".format(input_lyr_name)
    sys.exit()

## the mean center (average x and average y coordinate)
avg_x, avg_y = np.mean(xy_arr, axis=0)

print "Mean Center: {0}, {1}".format(avg_x, avg_y)

sum_of_sq_diff_x = 0.0
sum_of_sq_diff_y = 0.0

for x, y in xy_arr:
    # (x - xmean)squared
    diff_x = math.pow(x - avg_x, 2)
    # (y - ymean)squared
    diff_y = math.pow(y - avg_y, 2)
    # sum the differences sqaured from above
    sum_of_sq_diff_x += diff_x
    sum_of_sq_diff_y += diff_y

# x axis length
sum_of_results_x = (sum_of_sq_diff_x/lyr.GetFeatureCount())
standard_distance_x = math.sqrt(sum_of_results_x)
# y axis length
sum_of_results_y = (sum_of_sq_diff_y/lyr.GetFeatureCount())
standard_distance_y = math.sqrt(sum_of_results_y)

print standard_distance_x, standard_distance_y


texts: Snap existing line vertexes to nearest line in MapInfo

texts:
I have driven along some roads and registered some information about the state of the edge of the road. The GPS coordinates are not very precise, and I would like to smooth or snap the tracked lines to the actually edges of the road.
All data is registered driving at the right direction in the right side of the road. I have three layers with data.

Red layer is the edges of the road.
Blue layer is the center lines of the road.
Green layers is the tracks registered by the GPS with some attributes about the state of the side of the road. Look at the image at the bottom.

I have many of these and would like to automatically smooth or snap the points of the green lines to the right side of the road, to make them follow the red line.

How can I do that using MapInfo?

Just some mathematical description how to solve it would be nice. I am not asking for a full MapBasic script that solves it.

My initial thought is to find the nearest red line to a chosen green line. For each vertex in the green line get a vector to the nearest point on a red line segment and move the vertex by this vector. But I can see the GPS is sometimes very unstable and the green line is on the wrong side of the center line. In some way I need to find the nearest red line to the "right" of the blue center line and move the vertices to this.

I put some white arrows on the image to illustrate where to move the green lines to.

Perhaps a totally different way can be used?


texts:
Could you please explain the meaning of this equation:

(b1 le 0) * 0 + (b1 ge 10000) * 1 + (b1 gt 0 et b1 lt 10000) * float (b1) / 10000 b1.

What do ge, gl etc mean?


texts: How to create labeling similar to this map?

texts: Any ideas on if ArcMap is able to create labeling similar to this map? Ie. 95xx or 9500 is listed outside of the units and the the actual last 2 digits are listed within so that each unit does not have 4 digits labeling inside?


texts: I need help with merging two different tables in QGIS. First table is normal atribute table for spatial layer. Second table is text .dbf file without geometry. I cant use "join" feature, because in text file I have more (two or three) lines to one line in spatial layer. In ArcView 3 is "Link" feature and in ArcGis 3.2 is "Relate" feature. In QGIS I am helpless.


texts:
I don't know what's going on. My project CRS and layer CRS are fine, although if I export them to the KML file, I have something weird in my Google Earth.

I tried also to export to Shapefile and add this layer as a new one, but it didn't work either.


texts: Why is my PostGIS transformation not working?

texts:
I am using a PostGIS query to find the intersection of two streets, and then transform the results into lat/lon. The street database's current SRID is 6856, a projection I added to PostGIS myself; and the target SRID is 4326. For reasons unknown, the transformation is not working. Here is my original query, which produces a point in feet:

SELECT distinct(AsText(intersection(b.geom, a.geom))) AS the_intersection FROM (SELECT geom FROM public.streets WHERE streetname='LOVEJOY' AND prefix='NW' AND ftype='ST') a, (SELECT geom FROM public.streets WHERE streetname='23RD' AND prefix='NW' AND ftype='AVE') b WHERE a.geom && b.geom AND intersects (b.geom, a.geom);

And here is my modified query, in which I attempt to transform the feet to lat/lon:

SELECT distinct(AsText(ST_Transform(intersection(b.geom, a.geom),4326))) AS the_intersection FROM (SELECT geom FROM public.streets WHERE streetname='LOVEJOY' AND prefix='NW' AND ftype='ST') a, (SELECT geom FROM public.streets WHERE streetname='23RD' AND prefix='NW' AND ftype='AVE') b WHERE a.geom && b.geom AND intersects (b.geom, a.geom);

This is my error message:

ERROR:  AddToPROJ4SRSCache: could not parse proj4 string ''

I have verified that both 6856 and 4326 are available in my spatial_ref_sys. Any insights would be appreciated.


texts:
I'm trying to use polygons in arcpy. The polygons are stored in a shapefile. When running the following:

import arcpy
gdb = 'path//to//gbd'
inshape = os.path.join(folder,'Ameland_1968.shp')
poly = os.path.join(gdb,'test')

arcpy.MakeFeatureLayer_management (inshape, poly)

I get this error:

ExecuteError: ERROR 000229: Cannot open inshape **(red: this path is typed out)
Failed to execute (MakeFeatureLayer).

What is the problem? Arcpy should be able to run this line with shapefiles right? I'm actually trying to run a part of a script that's using arcpy but the data is in shapefiles. How to convert these?

edit: this is my original code

poly_shape = 'D:\\ARNO\\Analysis\\ArcGIS\\ArcPyScripts\\Data\\straight_bars_2\\Ameland\\Ameland_1965.shp'
poly = 'D:\\ARNO\\Analysis\\ArcGIS\\ArcPyScripts\\Data\\test.gdb\\Ameland65'
arcpy.MakeFeatureLayer_management (poly_shape, poly)


texts: How to correctly configure OSGeo4W64 for Windows 7 64-bit?

texts:
I have successfully installed OSGeo4W64 on numerous new systems (actually EC2 instances) and had no issues.  But I've run into an issue with it on a Windows 7 workstation that has  a lot of JREs on it.  Because of that, I set my paths from a batch to whatever JRE I need to use. But I can't get any of 64-bit versions to work:a ogrinfo --help throws the missing jvm.dll machine error message.

Which JRE should I be using and how should the pathing for the JAVA_HOME and PATH variables be set?

Update:

Removed all JREs/SDKs.
Installed JRE 1.7 (jre1.7.0_67) - 64 bit.
Set the JAVA_HOME and PATH variables.
Started ogrinfo.exe and received a different error:


texts: How is QGIS opening a shapefile with no *.prj file?

texts: I was wondering how QGIS could open a shapefile with no *.prj file (for instance, let's say we have a shapefile with just *.shp, *.dbf and *.shx files). If the On the fly reprojection option is activated, the layer is loaded and displayed successfully. I don't understand how such a reprojection is possible if no initial CRS is set.


texts: Results of volume calculation differ between cut/fill vs polygon volume in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.3 for Desktop.

We just have a project that require to create a pond. So I get to the field to collect survey points for topography and work expense purpose (we use volume of the land that had been removed to calculate the expense). I want to know the volume of land that were taken out so I created TINs and convert them to Raster (one is for the land before removed and another one is the raster of the land after it have been finished) and run cut/fill then the result give me about 57000 cubic metres but when I run Polygon volume (3d analyst) it gave me up to 82000 cubic metres for the same survey points and TIN interface. I used to use Polygon volume to calculate volume of reservoir of water but never experience with cut/fill but one of my friends told me that cut/fill method is more accurate with earthwork. As the results are not even near and I don't know the process of each function very well. I can't tell which one is the correct one. I am confused now.

Please help explain me with details.


texts: Visual Impact Assessment in QGIS or GRASS

texts: Is it possible using QGIS or GRASS to perform a Zone of Visual Impact assessment? E.g. where can a 100m high structure be seen at a viewing height of 1.8m within 30km, taking into account elevation (also maybe buildings and wooded areas)? I'm pretty new to GIS but am having to learn pretty quickly for work. Any help would be much appreciated.


texts: deleting Broken Source from arcgis toc with python addin

texts:
i have created one python addin tool. which create some layers in Memory. at the end of the script it deletes the In memory layer. but such layer still remains in TOC with a broken link

how can i remove such layer in python ?


texts: Vegetation Change Tracker and GEE: What bands are required?

texts:
I'm hoping to use the GEE implementation of the Vegetation Change Tracker (VCT) algorithm here: https://developers.google.com/earth-engine/apidocs/ee-algorithms-temporalsegmentation-vct#javascript
I'm struggling to identify the band requirements for the timeSeries argument. The docs indicate the following:

B3
B4
B5
B7
thermal
NDVI
DNBR
COMP

Bands 1-6 are straight forward enough but 7 and 8 are confusing. DNBR (7) is differenced Normalized Burn Ratio I think, but I would expect to pass NBR for a time series... What is COMP (8)?


texts: Is it possible to move one feature in front of another in the same layer?

texts: Is it possible to move feature in front of another? If I have smaller polygon and I draw larger polygon over it, is there a way how to move the smaller one in front? Is it possible in QGIS or better in PyQGis?


texts: Model Builder: Iterating files; Don't know what to make "input" in later tool

texts:
I have a simple model where I am trying to iterate through a year of daily tif's and recode everything in them as either snow covered ground or not snow.

My problem is I have no clue what to put into the "input raster" field when I get to the reclass by table tool. I'm limited on model builder and python experience and am walking myself through this using grad school tutorials I found. In some of these resources, its indicated to leave the input raster field blank as the iterate function will generate new names for each file it sends to the tool. This can't be the case since, as the photos show, the model isn't even complete without this input. I have also tried using the %Name% variable (model wouldn't verify) and also got desperate and just put something random hoping the iteration would somehow override it (stupid, I know).

What am I missing? Is there a larger issue in the way or is it as simple as just being new at this and not knowing the iterate output?

First image: dialogue for the iterate function showing some constraints on file names and types.
Second: dialogue for the reclass by table tool showing empty "input raster" field.


texts:
My question is similar to slow import via osm2pgsql to postgresql database and Optimizing osm2pgsql imports for OSM data but as we currently have quite some problems importing a large set of OSM data I open a new one.

What is the best way to import a large dataset (OSM Export of Europe) into a postgres DB?

Our computer has 32 GB of RAM... so it could use all of that.

We tried a couple of params, but had no success... last try we used the

osm2pgsql -c -S /usr/share/osm2pgsql/default.style --slim -d osm-europe -U postgres -C 25000 europe-latest.osm.pbf

But we ran out of memory even though our sever has 32 GB of RAM available.

pending_ways failed: out of memory for query result
(7)
Error occurred, cleaning up

How do we improve our import command?

Even if it takes longer... But we need to have the data imported into our postgres DB.

Would you recommend using an EC2 for the task or should our setup work with different parameters?


texts:
I received a series of *.img files with a single .hdr file that contains the necessary georeferencing information (i.e., samples, lines, coord system, etc.). I presume this is an ENVI header file since the first lines in the .hdr file are

ENVI
description = {
...

The img files display in ArcGIS with no problem, but I can't for the life of me figure out how to get Arc to recognize the header file, or how to set this information in the raster after its loaded. Ideally, it would be great if I could pull the entire series of img files into a single catalog.


texts: How to change the extent of a similar raster in QGIS?

texts:
I have two MODIS evapotranspiration rasters and their information is given below.
Desirable raster extents:

Erratic raster extents:

I want to change the erratic raster extent to desirable raster extent. How can I do this?


texts: I'm looking for options to add aerials to my QGIS Project. The OpenLayers plugin doesn't work (it says it's only available for versions 2.0 - 2.2), so I'm stuck moving from ArcMap to QGIS. Is there maybe a repository I could go to to download georeferenced aerials? I'd prefer not to have to go through Google Earth downloading a few hundred aerials and then georeferencing them. Or maybe there is a plugin I'm missing?


texts: Append data with ogr2ogr when target has two geometry columns

texts:
I am running the following command:

ogr2ogr.exe -update -append --config PG_USE_COPY YES -f PostgreSQL "PG:host=xxxx port=5432 dbname=yyyy password=***** user=zzz" "C:\xxx.gpkg" -where "Owner_ID in (89061,720883,713699,88665,8985008,708759,713696,721202,721203,715717,9394015,9394258,9393828,720197,719592,721195,718608,9379306)" -lco LAUNDER=YES -lco GEOMETRY_NAME=geom -nln schema.existing_table -a_srs EPSG:27700

The attribute columns in the table have the same names (albeit in caps) but the PG table has more columns

All the attributes are copied across correctly, but I can't work out how to specify where the geometry should go as there are two geometry columns in the PostGIS table. I don't think any of the layer creation options work because it's not actually creating a layer.

Is this currently possible?


texts: Encoding properties for each point in a GeoJSON multipoint

texts:
I'm working with GeoJSON Multipoints and wondering if there's any way of encoding properties for each point in the array? I realise this is a vague question but I'm looking for suggestions as to how I might achieve this or if there's a standard I'm unaware of.

AN EXAMPLE

In the example below I would be looking for a way for point1_property to represent a value for the first point in the MultiPoint array and point2_property to represent a value for the second point.

{
    "type":"FeatureCollection",
    "features":[
        {
            "type":"Feature",
            "geometry":
                {
                    "type": "MultiPoint",
                    "coordinates": [
                        [100.0, 0.0], [101.0, 1.0]
                    ]
                },
                "properties":{
                "point1_property":"value1",
                "point2_property":"value2"
            }
        }
    ]
}


texts: Translate postgis data to a coordinate system

texts: I have a PostGIS database where the data is in an unknown or incorrect coordinate system.  Is there a way given some known points in the data to translate all the data to a specific coordinate system.


texts: Is there a smartphone app for mobile asset logging?

texts:
I'm after a recommendation for a smartphone app (iPhone or Android) to log asset information for litter bins. I need the following features

GPS
geotagged photos
Form with drop down menus for user definable attributes, e.g. construction material,classification, condition etc
export to csv or shp


texts: How to create multi-ring buffers automatically in QGIS?

texts: I'm wondering if anyone knows of a tool or plugin for QGIS that allows for multi-ring buffers.  I've found a post about doing multiple buffers one by one, but was hoping there was a way to do this all in one place.  I know in ArcMap that you can set the number of buffers you want and distance, so just hoping there's something out there in QGIS to accomplish this.


texts: QGIS renders enclosed polygon wrong compared to ArcGIS

texts:
I have a street pavement polygon shapefile that is enclosed and hollow on the inside. In QGIS the entire polygon becomes a single blob. The file was downloaded from a county FTP site and works fine in ArcGIS but never appears correctly in QGIS. I don't know how to correct the geometry.

The first image is how it correctly appears in ArcGIS and the second is how it renders in QGIS.


texts: Clip image by each feature in a FeatureCollection with Google Earth Engine

texts:
I have a raster image myRaster (ee.Image) and a FeatureCollection indexGrid which contains a polygon grid with a single attribute index (string).

My goal is to split the image by polygons in a FeatureCollection and export each tile to Drive. The best thing I came up with follows, but it doesn't work:

var grid_export = function(feature){
  Export.image.toDrive({
  image: myRaster,
  description: 'part_' + feature.get('index'),
  folder: 'dst',
  scale: 30,
  region: feature.geometry()
  });
  return null;
};

indexGrid.map(grid_export);


texts:
I am trying to calculate the coverage that a set of points has of a road network layer. I have 100 million points and approximately 80,000 segments in my road network layer.

I have created a buffer around each road segment (roadsegment.shp) and am attempting to run a spatial query (using the standard spatial query tool in QGIS) to see which segments contain at least one of my 100 million points (points.shp). However, this query is taking far too long to run and I need to find a more efficient way to do this.

I am not concerned with the number of points in each road segment buffer and I do not need to pass an attribute from the buffer to ALL of the points. I just need to know whether or not a segment contains a point. Once I know that it does, it can be dropped from the analysis. I think that this must help me to make the query more efficient.

I should point out that I am fairly new to QGIS and only know how to use the standard tools in the software.

Can anybody help me with a script to speed up the processing?


texts: Why the outcome of $area (geometry), is different from area saga shape indices tool in QGIS?

texts:
I am puzzled on why I get different outcomes on trying to calculate the area (m2) for multiple polygons in QGIS.

My polygons consist of tiny objects of few cm extent.
When using $area (field calculator / geometry) for the smallest polygons I get 0.0000 while saga tool gives a value.

It seems saga shape indices area calculation is more robust, but is there a reason why $area (geometry) fails to produce an outcome and even in visually different polygons it gives similar values or even bigger values in smaller polygons?

The EPSG is set to 32631 31N (meters) so it shouldn't cause any trouble.
I don't know what I am missing here.

$Area (as calculated by $area), area (add geometry column Project CRS), Area (saga shape indices tool)


texts: Generating shapefile containing UTM grid (graticule)?

texts:
I'm looking for a shapefile containing an 1000m or 2000m UTM grid ("graticule"). All I was about to find on the web is a 10000m UTM grid (here: http://en.giswiki.net/wiki/UTM) and I am aware that a 1000m UTM grid for the whole world might be quite big.

That's why I'm now wondering if and how it would be possible to generate such a shapefile for a part of the world. I would prefer to use open source tools like gdal. Others output formats than a shapefile might be ok too.

Has anyone an idea?


texts: Publishing Geoprocessing Service with importing of numpy in ArcGIS Server?

texts:
I use a lot of scripts as GP Service almost all use the arcpy functions but one (and very important) of my scripts needs a gdal/ogr and numpy imports. I can run this script in PyScripter or as a tool in Toolbox in ArcMap. Now, I am trying to publish my script from GP Result window as GP Service on ArcGIS Server 10.1. This is weird because I can see only Packaging in progress window a long time and after 30 minutes I must kill the Service Editor.
I tested my script by commenting it line by line and publishing each time as a GP service. I have already found the reason of my problems but still I do not have solution.  The publishing problems appear when I use import of numpy  in my script (import numpy as np or import numpy). Does anyone know the reason of this state? Numpy is very important library so I think the solution should exists.

I have tested my simple script

import arcpy
from arcpy import env

from osgeo import gdal
from osgeo import ogr
from osgeo.gdalconst import *
import numpy as np

arcpy.AddMessage("End")
print("End")

It runs in ArcGIS for Desktop but I can not publish service. I always get window Packing in progress

When I comment out import as np I can publish service successfully.


texts: GeoJSON is valid according to JSONLint but jquery gives error

texts:
I have the following GeoJSON (from webservice)

{
"Features": [
    {
        "Id": "1",
        "Geometry": {
            "Coordinates": [
                {
                    "Latitude": 31.2812,
                    "Longitude": 34.8027,
                    "Altitude": 0
                }
            ],
            "Type": 0,
            "CRS": null,
            "BoundingBoxes": null
        },
        "Properties": {
            "Color": "#00FF55"
        },
        "Type": 7,
        "CRS": null,
        "BoundingBoxes": null
    },
    {
        "Id": "2",
        "Geometry": {
            "Coordinates": [
                {
                    "Latitude": 31.7825,
                    "Longitude": 35.1966,
                    "Altitude": 0
                }
            ],
            "Type": 0,
            "CRS": null,
            "BoundingBoxes": null
        },
        "Properties": {
            "Color": "#FF5500"
        },
        "Type": 7,
        "CRS": null,
        "BoundingBoxes": null
    }
],
"Type": 8,
"CRS": null,
"BoundingBoxes": null

}

This is the webservice code :

public string GetItemsFromDBGeoJSON()
{
    List<GeoJSON.Net.Feature.Feature> fList = new List<GeoJSON.Net.Feature.Feature>();
    using (SqlConnection conn = new SqlConnection(connectionString))
    {
        using (SqlCommand cmd = new SqlCommand("dbo.GetGPSData", conn))
        {
            conn.Open();
            SqlDataReader reader = cmd.ExecuteReader();

            while (reader.Read())
            {
                string name = (string)reader["GPSUserName"];
                double lon = Convert.ToDouble(((decimal)reader["Lon"]));
                double lat = Convert.ToDouble(((decimal)reader["Lat"]));
                string clr = (string)reader["GPSUserColor"];

                GeoJSON.Net.Geometry.GeographicPosition latlon = new GeoJSON.Net.Geometry.GeographicPosition(lat,lon,0);
                GeoJSON.Net.Geometry.Point point = new GeoJSON.Net.Geometry.Point(latlon);

                Dictionary<string,object> properties = new Dictionary<string,object>();
                properties.Add("Color", clr);
                GeoJSON.Net.Feature.Feature feature = new GeoJSON.Net.Feature.Feature(point, properties);
                feature.Id = name;

                fList.Add(feature);
            }
            conn.Close();
        }
    }
    GeoJSON.Net.Feature.FeatureCollection fCol = new GeoJSON.Net.Feature.FeatureCollection(fList);
    JavaScriptSerializer serializer = new JavaScriptSerializer();

    return serializer.Serialize(fCol);
}

The Ajax call gives me error :

SyntaxError: JSON.parse: unexpected character

Here is the request using JQuery

$.ajax({
    type: "POST",
    url: "http://localhost:36361/Service.asmx/GetItemsFromDBGeoJSON",
    dataType: 'json',
    success: function (response) {
        console.log('success');
    },
    error: function (req, status, error) {
        console.log(error);
    }
});

Also , Using Firebug Console , I see this response from the request:

<?xml version="1.0" encoding="utf-8"?><string xmlns="http://tempuri.org/">{"Features":[{"Id":"1","Geometry":{"Coordinates":[{"Latitude":31.2812,"Longitude":34.8027,"Altitude":0}],"Type":0,"CRS":null,"BoundingBoxes":null},"Properties":{"Color":"#00FF55"},"Type":7,"CRS":null,"BoundingBoxes":null},{"Id":"2","Geometry":{"Coordinates":[{"Latitude":31.7825,"Longitude":35.1966,"Altitude":0}],"Type":0,"CRS":null,"BoundingBoxes":null},"Properties":{"Color":"#FF5500"},"Type":7,"CRS":null,"BoundingBoxes":null}],"Type":8,"CRS":null,"BoundingBoxes":null}</string>

Any idea why ? I want to show it on leaflet layer.


texts: Adding WMTS layer to Mapbox GL JS from GeoServer

texts:
I'm trying to add WMTS layer from GeoServer 2.18 in Mapbox GL JS. My code:
            window.mymap.addSource(
            "mvt-test-source", {
             "type": "vector",
            "tiles": [
                "http://localhost:8080/geoserver/gwc/service/wmts?SERVICE=WMTS&&VERSION=1.0.0&REQUEST=GetTile&layer=mapbox:roads4326&TILEMATRIX=EPSG:4326:{z}&TILEMATRIXSET=EPSG:4326&FORMAT=application/vnd.mapbox-vector-tile&TILECOL={x}&TILEROW={y}"
                ],
            "minZoom": 0,
            "maxZoom": 14
        });

        window.mymap.addLayer(
            {
                'id': 'my_mvt_layer',
                'type': 'line',
                'source': 'mvt-test-source',
                'source-layer': 'mapbox:roads4326',
                "visibility": "visible",
                'paint': {
                    'line-color': 'Red',
                    'line-width': 7
                }
            }
        );

There are no errors in browser, but layer does not appear on the map. What could be the reason?
May be Mapbox GL JS does not support this protocol? When I try to use WMS service - all work fine. Next code with WMS format:
window.mymap.addSource("wms-test-source", {
                "type": 'raster',

            "tiles": [
              "http://localhost:8080/geoserver/wms?bbox={bbox-epsg-3857}&format=image/png&service=WMS&&version=1.1.1&request=GetMap&srs=EPSG:900913&transparent=true&width=256&height=256&layers=mapbox:roads4326"
            ],
            'tileSize': 256
        });

        window.mymap.addLayer(
            {
                "id": 'wms-test-layer',
                "type": 'raster',
                "source": 'wms-test-source',
                "paint": {
                }
            }
        );


texts:
I'm working on a project in India and have problem with match shapefile to Google base map. I have tried different projections but couldn't find the right one.
What is the coordinate system to use in Ramthal area?


texts: Dissolving streets by distance using ArcGIS Desktop

texts:
I am using ArcGIS Desktop.
Is there any way to dissolve a streets network to districts based on a distance given?
The streets of the districts must have the same sum of the  distance.
example:
The selected districts should look like the one in blue, the sum of the streets in this district is 2 km


texts: Creating parallel polylines with a constant distance

texts: Can ArcMap fill in selected polygon with many parallel polylines with a constant distance between them and with a particular azimuth?


texts: Spatial join Polygons by their centroid in ArcGIS

texts:
I am trying to spatially join two tables containing polygons (let's call them table A and table B) so the ID of table A is joined to the table B only if the polygon centroids of table A are inside the polygons of table B.

Hope it makes sense?

I tried the spatial join, but the only option I could select is if the polygons of A are completely within the polygons B (second radio button) which is something I don't want.


texts: Openlayers 3 ol.layer.Tile shift from version 3.6.0 and forward

texts:
I have local tiles on my OL 3 app using version 3.6.0.
I've noticed when I upgrade to version 3.7.0 I need to rewrite my tileUrlFunction or else they don't shown anymore.

This is the code change I need to do :

 source: new ol.source.XYZ({
    tileUrlFunction: function (coordinate) {
        if (coordinate == null) {
            return "";
        }
        var z = coordinate[0];
        var x = coordinate[1];
        //var y = coordinate[2]; //3.6.0
        var y = -coordinate[2]; //3.7.0

        return 'tiles/controlroom/' + z + '/' + x + '/' + y + '.png';
    },
    minZoom: 7,
    maxZoom: 16
})

Once I change to 3.7.0 there is big deviation the further my zoom out is

Example Max Zoom Out:

When I zoom In , the gap seems to get closer where it used to be with 3.6.0

Why does this happens ?


texts: mongodb $geoIntersects query for all documents who's polygon contains a point. no results?

texts:
I am trying to  find all documents for a point that is within a document's polygon.

using mongoTemplate to set up my geo index.

 mongoTemplate.indexOps(Boundary.class).ensureIndex(
        new GeospatialIndex("polygon"));

I have a Criteria2ds that extends Criteria:

public class Criteria2ds
 extends Criteria {

    private Map<String, Object> criteria = new LinkedHashMap<String, Object>();
    ...
    ...
    public static Criteria2ds where(String key) {
        return new Criteria2ds(key);
    }
    public Criteria2ds intersects(Point point) {
        Assert.notNull(point);
        final ArrayList<Double> coordinates = new ArrayList<Double>();
        coordinates.add(point.getY());
        coordinates.add(point.getX());
        BasicDBObject geometry = new BasicDBObject()
            .append("type", "Point")
            .append("coordinates", coordinates);

        criteria.put("$geoIntersects", new BasicDBObject("$geometry", geometry));

        return this;
    }

my query:

Query   query = new Query();

query.addCriteria(Criteria2ds.where("polygon").intersects(point));

..
this creates the query:

    { "polygon" : { "$geoIntersects" : { "$geometry" : { "type" : "Point" , "coordinates" : [ -105.01621 , 39.57422]}}}}

Collection<MyDoc> list  =  mongoTemplate.find(query, Boundary.class)

however, the results are 0?

is there something wrong with my code?

public class MyDoc
       Polygon polygon;
       public Collection<Polygon> getPolygons() {
        return polygons;
    }

    public void setPolygons(Collection<Polygon> polygons) {
        this.polygons = polygons;
    }

public interface BoundaryRepository extends MongoRepository<Boundary, String>{
}

I save in a DAO like:

boundaryRepository.save(myDoc);


texts: "slyr-master" plugin for QGIS doesn't work due to "missing meta data"

texts:
I downloaded the "slyr" plugin from GitHub which should be able to convert .style data into QGIS. According to the README file, I'm just supposed to add the "slyr" folder to the Python -> Plugin folder and activated it in QGIS.

When I go to the QGIS Plugin extentions the plugin appears but I get the error message "This plugin is defect. Metadata is missing". However, there is a metadata.txt file inside the "slyr" folder.

Does anyone have any ideas what the problem might be and how to solve it?

I'm using QGIS 3.4.6.


texts: Can raw gps data for each individual satellite be extracted from android phones?

texts: I am trying to build a tight integration module between GNSS and INS. Hence i need raw satellite data. Can this be extracted from normal android phones?


texts: Postgis performance increase with cluster?

texts:
I've been reading a few questions on this site and saw quite a few people advocating for creating a cluster function.

So I've decided to test this function to see if there was any notable speed increase. I've used a table with 35,000 points. Both tables have the spatial index but tableb has the cluster

create table tableb as select * from tablea

CREATE INDEX tableb_sidx ON tableb  USING gist  (geom);

CLUSTER tableb USING tableb_sidx

select a.*
from tablea a, inventory b
where st_dwithin(a.geom, b.geom, 200)

select a.*
from tableb a, inventory b
where st_dwithin(a.geom, b.geom, 200)

Running these st_dwithin queries result in the exact same running time of 1 minute

Am I doing something wrong or does this cluster not apply to every spatial query?


texts: Arcobjects .Net, suppress geoprocessing messages

texts:
I am executing the copy features tool via the geoprocessor inside my .net custom geoprocessing function tool. When the copy features tool executes it is writing it's messages to my tool's results. I want to stop these messages from being written and handle them myself. Is there a way to stop a geoprocessor from writing messages? I am writing .net on 10.1.

This is the code I am running the tool with.

            IGeoProcessor2 gp = new GeoProcessorClass();
            IGeoProcessorResult results = new GeoProcessorResultClass();
            IVariantArray parameters = new VarArrayClass();
            parameters.Add(tempFC);
            parameters.Add(outFC);
            results = gp.Execute("CopyFeatures_Management", parameters, null);


texts: Is there a way to create a plugin in C++ without compiling all the Qgis source?

texts:
This is my second post and I apologize if it's post the same day on the first.

After several researchs, I found the only way to create a plugin in C++ is to compile the Qgis source with the plugin created thanks to the "plugin_builder.py" which can be found in the "src/plugins" on the git qgis (https://github.com/qgis/Quantum-GIS/tree/master/src/plugins)

To create a simple plugin which can be add simply by putting the folder of the plugin
in ".qgis/python/plugins" or "C:/OSGeo4W/apps/qgis/python/plugins" the plugin must be written in Python, It is rigth ?

So my question is simple :
Is there a way to create a plugin in C++ without compiling all the Qgis source ?

Thank you very much
Regards


texts: Layout Composer Text Box Expression - isolating names from project file path

texts:
In QGIS 3.X, I'm attempting to use an expression in layout composer to automatically input the client and project name into a text box. The details are always present in the file path to the project as you can see below. The number of backslashes will not change and the client name will always be between the 8th & 9th. Likewise, the project name follows the project number (which should be omitted from the expression.) A snip from the layout composer of how the text box should look is below.

C:\Users\user\Dropbox\Encrypted\GIS\1. Client Projects\S\Smith, S\ABC999999 - Site Finding, England


texts: Are there any geocoding plugins for QGIS that interpolate along linear features instead of querying web services?

texts:
I am curious if there are any tools for QGIS to facilitate geocoding tabular data using a local linear features layer (ie, a street map with street name and address ranges - similar to geocoding in ArcGIS). I am aware of Alessandro Pasotti's GeoCoding plugin, but it queries Google Maps/OpenStreetMap rather than interpolating along line segments of a user-provided layer.

One reason I am interested in this approach is that for some projects it is desirable to have geocoded points placed exactly on features from other project-specific layers, such as the street maps provided by local services or municipalities, which may not necessarily match the positions returned by web services.

Any ideas - or any input in case I have a go at writing a solution from scratch?


texts: QGIS distance from a cell center to the nearest point

texts:
I am trying to run the following analysis in QGIS 3.2.

My inputs:

raster layer (with population estimates from LandScan)
vector layer, with points (each point represent the location of a mobile cell tower)

I want to define for every cell of the raster layer (masked for one specific country):

it's value
it's distance (preferably in meters) to the closest point/cell tower.

The output can be another raster layer or any format that can be exported for further analysis in Excel.

I have tried so far Raster -> Nearest Neighbour, which generates a binary raster map with no values, so I am clearly doing something wrong.

A snapshot of my current setup is below. Blue dots are cell towers, a brown gradient is the raster data.


texts: Determining if shapefile and raster overlap in Python using OGR/GDAL?

texts:
I am building a script in python using OGR/GDAL.

I have a set of shapefiles and a set of GeoTiff raster files.

I would like to have my script ignore shapefiles if they do not intersect with the raster area.

The shapefile is not a rectangle, so I cannot simply compare the xmin/xmax,ymin/ymax values returned by layer.GetExtent(). I need the actual polygon representing it's overall shape, and then some way of determining if that polygon intersects with the raster square.

I was thinking I could somehow merge all the polygons in the shapefile into one feature, and then read the geometry on that feature, and then compare that information to the raster extent. However, I am unsure of specifically how to execute this.

How to extract border polygon information from shapefile?
How to determine if that polygon intersects a given square area?


texts: How can I modify the tiles size in Openscales or other ways to boost the application performance

texts:
I am developing an application in Flex that uses Openscales and GeoServer in order to provide a map in the background. It works but the perfs aren't really good; even if everything is on the same computer (the application, the GeoServer and the MySQL database or the shapefiles that are used by GeoServer), loading the map takes some time (~5-10s for the whole screen).

I thought that with everything local, it would be much more fluid than that. For example, when you drag the map, I'd like the map to be updated almost instantly.

In order to boost this, I tried many things on Geoserver (shapefiles, MySQL, ...) and modified many parameters but it doesn't change anything. I'm now trying to change the tiles size : I think that they are for the moment in 128x128 so I would like to change their height and width so that the app asks the server for less tiles each time the user drag the map.

My questions are:

How do you change those width and height? I coudn't find a way to do it.
How could I increase my app performance regarding the map?

The code that creates the map object is the following :

<?xml version="1.0" encoding="utf-8"?>
<s:Group xmlns:fx="http://ns.adobe.com/mxml/2009"
     xmlns:s="library://ns.adobe.com/flex/spark"
     xmlns:mx="library://ns.adobe.com/flex/mx"
     xmlns:os="http://openscales.org"
     xmlns:layer="org.openscales.fx.layer.*" xmlns:core="org.openscales.core.*">
    <os:Map id="fxmap" width="100%" height="100%" resolution="0.01" maxExtent="-180,-90,180,90" projection="EPSG:4326" center="0,45">
        <os:WMSC id="mapLayers"
             name="Topp States (WMS)"
             url="http://localhost:8080/geoserver/gwc/service/wms"
             layers="WorldOnMySQL:ne_10m_admin_0_countries"
             format="image/png"
             bgcolor="0x454E54"
             projection="EPSG:4326"
             tileHeight="256"
             tileWidth="256"
             tiled="true"
             />
        <os:DragHandler/>
        <os:ClickHandler/>
        <os:WheelHandler/>
    </os:Map>

After that, there isn't anything useful in the script part...

I hope you can help me, and ask me if you need more information.


texts: How to make a SQL query that would sum multiple records in ArcGIS?

texts:
I've got a large file geodatabase that contains a polyline feature class of water-pipes. The feature class have many fields, but the ones that I'm interested in is that of the water-pipe's Length and Locality Name.

I'm looking for a way to use the Table Select tool to generate a summary table in the same file geodatabase which summarizes the water-pipe's shape lengths for each locality. The output should be file geodatabase table which looks like this example:

Locality 1      24395.66118
Locality 2      19572.35849
Locality 3      17474.87241
Locality 4      16205.04598
etc...

I've tried to follow this link which describes the SQL syntax for such operation. I didn't get where to write it and how. This is the syntax I probably need from the link:

SELECT STATE_NAME, SUM(POP1990) as TotalPopulation
FROM counties
GROUP BY STATE_NAME ORDER BY STATE_NAME


texts: UTM Projection of DEM-data with gdalwarp

texts:
I have downloaded a certain tile from https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm, and after the requisite logistics (unzipping, putting all *.AVE_DSM.tif files into a single folder), I built a virtual
dataset with gdalbuildvrt *.tif asia.vrt.
Now gdalinfo asia.vrt gives this:

Size is 108000, 90000
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (89.999999999999957,50.000000000000028)
Pixel Size = (0.000277777777778,-0.000277777777778)
Corner Coordinates:
Upper Left  (  90.0000000,  50.0000000) ( 90d 0' 0.00"E, 50d 0' 0.00"N)
Lower Left  (  90.0000000,  25.0000000) ( 90d 0' 0.00"E, 25d 0' 0.00"N)
Upper Right ( 120.0000000,  50.0000000) (120d 0' 0.00"E, 50d 0' 0.00"N)
Lower Right ( 120.0000000,  25.0000000) (120d 0' 0.00"E, 25d 0' 0.00"N)
Center      ( 105.0000000,  37.5000000) (105d 0' 0.00"E, 37d30' 0.00"N)
Band 1 Block=128x128 Type=Int16, ColorInterp=Gray

I can query the data as follows:
gdallocationinfo -geoloc asia.vrt 110.00 30.00, meaning "Find the elevation
of the point at 110 degrees East and 30 degrees North",
and get a proper

Report:
  Location: (72000P,72000L)
  Band 1:
    <LocationInfo><File>N029E110_AVE_DSM.tif</File></LocationInfo>
    Value: 894

Now I would like to project the entire region I've got here
to a rectangle P = [a,b]x[c,d], so that I can query not with lat/lon pairs but
rather with a point from P. Naturally, I need the distance
between the two points be somewhat reflective of their "original" distance.
I do not care much about the specific values of a,b,c,d, and anyway I think
they depend on the units option we would supply to gdalwarp.
Now my question is,
what is the exact gdalwarp-command I need to issue to accomplish this?
I know it has to do with proj-string, but I am lost to its multitude of parameters.


texts:
On Google Earth Engine I'm trying to download some display images using the method image.getThumbURL. The image shows well in Map.addLayer, but when I pass it through getThumbURL I only get a blank image. Here an example of this behavior:
var point = ee.Geometry.Point([11.300422012780208, 44.50570009326473])
var pngArea =
    ee.Geometry.Polygon(
        [[[10.966712295983333, 44.65340259358304],
          [10.966712295983333, 44.24754759759594],
          [11.738501846764583, 44.24754759759594],
          [11.738501846764583, 44.65340259358304]]], null, false)

var s2coll = ee.ImageCollection("COPERNICUS/S2")
s2coll = s2coll.filterBounds(point)
        .filterDate('2018-07-15', '2018-07-30')
print(s2coll)
var img = s2coll.mosaic()
print(img)
var visualization = {
  min: 0.0,
  max: 2000,
  bands: ['B4', 'B3', 'B2'],
}
var outImg = img.visualize(visualization)
Map.addLayer(outImg)

print(outImg.getThumbURL({region: pngArea, format: 'png'}))

https://code.earthengine.google.com/89afe6e672f63d2245ef438288f78f55
What am I doing wrong?


texts: Pyqgis: How to use GDAL-Tool 'Rasterinfo'

texts:
I have a couple of raster layers and for my next Plugin, I want to get the corner coordinates for each layer.

My idea was to use the GDAL-Tool 'Rasterinformation' which gives as an output a lot of information about a raster file including the corner coordinates. However, after running the tool, the output is presented in the GUI.

How do I manage this in PyQGIS? My code so far is the following:

    def selectRasters(self):
    inputDir = self.getInputPath()
    for lyr in glob.glob(os.path.join(inputDir, "*xyz")):
        info = []
        processing.run('gdal:rasterinfo', lyr, False, True, True, info)

My hope is that I have all the information as a list and then I can do something like raster upper left = info[5] and so on.

Long story short I want add the corner coordinates into a variable. Any ideas how to do that?


texts: Finding tight corners on roads using ArcGIS Desktop?

texts:
Problem:

How can I pinpoint tight corners on a road dataset that might cause longer/larger vehicles difficulty in driving around.

Assumptions:

Ignore turns into junctions.  Just interested in corners along a single road.
Road network is a polyline layer.  No arcs involved.
Road width is an attribute.
A large vehicle can be defined as a semi-truck.

Ideal Solution:

The end goal would be for a map showing potential areas on a road network that would be best avoided by larger vehicles.  Ideally this would have a score value as opposed to a yes/no value, to signify how tricky the particular corner is, depending on the size of vehicle.

I have likely not provided enough information, so feel free to ask more Qs via comments and ill supplement my Q with more info.  For info, this was a question I was posed, and I had some ideas, but would love to hear what your thoughts are on a potential workflow to solve this question using ArcGIS Desktop.


texts: NetCDF layer (.nc) appears flipped and rotated in QGIS

texts:
I'm trying to import NetCDF layers with environmental information to QGIS and they appear flipped and rotated

In the top-left you can see the Europe coastal line I was using, in the bottom-right the nc layer

When importing them to a friend's ArcGis I can choose which subdataset I want to use as X and Y coordinates, there are four candidates: lat, lon, nlat and nlon. After some trials and errors I finally got it right in ArcGis, but in QGIS importing it (Layer > Add layer > Add raster layer) there isn't any option to choose the coordinates.

After doing some research I found in a another post here a way of fixing the y-axis with GDAL in OSGeo4W

    gdal_translate -of netCDF -co WRITE_BOTTOMUP=NO NETCDF:"filename.nc":datasetname newfile.nc

For example:

gdal_translate -of netCDF -co WRITE_BOTTOMUP=NO NETCDF:"ATS_LST_3PUUOL200208XX_0.50X0.50_CS3_LSM.nc":LST 200208_05_LST.nc

Althought I got some warnings:

Warning 1: dimension #2 (nlat) is not a Longitude/X dimension.
Warning 1: dimension #1 (nlon) is not a Latitude/Y dimension.
Warning 1: dimension #0 (diurnal) is not a Time or Vertical dimension.
Input file size is 360, 720
0...10Warning 1: dimension #2 (nlat) is not a Longitude/X dimension.
Warning 1: dimension #1 (nlon) is not a Latitude/Y dimension.
Warning 1: dimension #0 (diurnal) is not a Time or Vertical dimension.
Warning 1: dimension #2 (nlat) is not a Longitude/X dimension.
Warning 1: dimension #1 (nlon) is not a Latitude/Y dimension.
Warning 1: dimension #0 (diurnal) is not a Time or Vertical dimension.
Warning 1: creating geographic file without lon/lat values!

Now the y-axis seems fine

Is there a way to do this for the x-axis? I found nothing here

This layer in specific is called ATS_LST_3PUUOL200208XX_0.50X0.50_CS3_LSM from CEDA

Here are the files in my dropbox if you don't have access to the webpage


texts:
I'm doing some tests with the v.clean but getting some results that I don't really know if I am the one doing something wrong or what.
As you can see in the image bellow I created 3 features, one with an area quite small compared to the others. I was trying to verify if the described here (https://www.oreilly.com/library/view/learn-qgis-/9781788997423/01a0973e-ce1d-474f-bcdc-2d4605c7b67a.xhtml) would work, and the feature with values bellow the threshold would merge with the neighboring polygon with the longest common boundary.
It simply did not happen. Why is that? Am I doing something wrong?
The area of the smaller polygon is 37.69 and the threshold set in v.clean was 40 with the cleaning tool rmarea selected.


texts:
I'm a newbie in gdal2tiles & python but I would like to modify gdal2tiles to create 512px tiles for retina display (screen).

Can someone suggest how I can do this?


texts:
I tried to use gdal_calc.py to avoid using the manual GUI raster calculator in QGIS. However, I'm not yet familiar with it and having errors with my code:
gdal_calc.py -A perennials.tif -B annuals.tif --outfile=annuals_only.tif --calc="(A<=85 or B=2001)*(A)" --format=GTiff 0 .. evaluation of calculation (A<=85 or B=2001)*(A)

evaluation of calculation (A<=85 or B=2001)*(A) failed Traceback (most recent call last):
  File "/usr/bin/gdal_calc.py", line 401, in <module>
    main()
  File "/usr/bin/gdal_calc.py", line 394, in main
    doit(opts, args)
  File "/usr/bin/gdal_calc.py", line 288, in doit
    myResult = eval(opts.calc, global_namespace, local_namespace)
  File "<string>", line 1
    (A<=85 or B=2001)*(A)
               ^
SyntaxError: invalid syntax

Is there some thing I can improve in my --calc formula?


texts: Creating mesh that divides world by same area?

texts:
How can I create a mesh that divides the world by the same area?

If possible, I would like to create it using qgis and postgres
I'm planning to use quadkey or geohash but If you already have standard
I want to use it


texts: GRASS iterate over multiple shapefiles using g.region

texts:
I have a shapefile called Ahmednagar. I have daily rainfall data from 1997 till 2008. I want to sum total rainfall from 152 till 273 day for each year for 1997 till 2008. This is how I do it.

g.region vect=Ahmednagar  #first set g.region as Ahmednagar
j=1997
while [ $j != 2008 ]
do
    r.mapcalc "tmpsum=0.0"
    i=152                   #start from 152 day of the year
    while [ $i != 274 ]     # till the 273 day of the year
    do
            g.copy rast=prec$j.$i,tmp #copies rainfall of jth year and ith day
            r.mapcalc "tmpsum2 = tmpsum+ tmp"
            r.mapcalc "tmpsum=tmpsum2"
            g.remove rast=tmp
            i=`expr $i + 1`
    done
    g.rename rast=tmpsum,precAhmednagar$j   #rename the tmpsum as the shapefile
    j=`expr $j + 1`
done

This will give me a raster files for each year from 1997 till 2007 called precAhmednagar$j where j is the year.

Now I have many more shapefiles (215 exactly). For each of them I want to do the same thing. But I am not sure how to go about it. Basically set a shapefile using g.region, do all the calculation and most importantly rename the final raster tmpsum as the precXXXX$j where XXXX is the name of the shapefile.

I have imported all the 215 shapefiles in GRASS the list which I can get using g.list vect but that's it.


texts: Finding regularly intervalled GPS positions between two known GPS points

texts:
I have been trawling for a solution for a while now, and either I cant find the correct solution or I'm not quite getting my search right.
I have some GPS points from a ship track that are collected every 10 seconds or so with the time they were collected. At one point the GPS goes down. I want to recreate estimated points during this down time. Is there a way to create a list of points with GPS coordinates between the last GPS fix and the next fix every 10 seconds?


texts: Debugging SQL in arcpy.da.SearchCursor RuntimeError: An invalid SQL statement was used?

texts:
I am going through a python script in arcpy 10.1 to Select a field value (i.e. NEAR_DIST) based on three condition as

        1) NEAR_DIST<=radius(i.e.500 or anything from user input)
        2) NEAR_DIST<>0
        3) NEAR_DIST<> -1

I am using 500 which is to be variable asked from user and to be saved in "radius" variable so, i am only using "radius" variable in the SQL statement in following code-

    import arcpy,os,sys
    feature =     "C:\Users\Bakul\Desktop\Od\example_shapefile\Polyline_segments_example.shp"
    radius = 500
    # create a dictionary to hold the list of near ids for each input id
    nearest_dict = dict()
    arra = []
    whereClause = '''"NEAR_DIST" <= '%{x}%' AND "NEAR_DIST" <>'%{y}%' AND "NEAR_DIST"   <>'%{z}%\''''.format(x=radius,y=-1,z=0)
    print radius
    with arcpy.da.SearchCursor(feature, "NEAR_DIST", whereClause) as rows: ##NEAR_DIST  is a field, I need sort all NEAR_DIST <=500 (i.e radius) and not equal 0 and -1
        for row in rows:## debugger does not pass here,
            nearest_id = row[0]
            arra.append(nearest_id)
        del rows
        del row
    print arra

Debugger Error is

    500
    Traceback (most recent call last):
      File "C:\Users\Bakul\Desktop\Od\Python\scratch.py", line 89, in <module>
        for row in rows:
    RuntimeError: An invalid SQL statement was used.

Can you tell where is the problem (i think in SQL whereClause)?


texts: How to track a brand new path with realtime NMEA data in Google Earth?

texts:
I have some trails on my property that are not visible on Google Earth or Google Maps.  I'd like to document them in Google Maps's "My Maps".  I believe I need to create KML files of the trails.
So my thinking is I connect my bluetooth GPS receiver to Google Earth.  I can do this with GE's realtime function and it shows my location.  I'll then walk my trails.  But not sure how to tell GE what I'm trying to do so that it outputs KML representing my trails.


texts: Setting basemap layers of WebMap using ArcGIS Python API?

texts:
Using the ArcGIS Python API, I've create a WebMap and two Layers (one containing features, the other containing labels):

cbmt = WebMap()
cbmt_geometry_layer = Layer(url = <url of geometry layer>)
cbmt_text_layer = Layer(url = <url of text layer>)

I'd like to add both of these layers to the map as basemap layers, like you can in ArcGIS Online or Portal when adding layers to a web map. The WebMap.basemap property is read-only. Is there a property on the WebMap.add_layer() method that specifies to add layers as basemap layers?

xref: https://community.esri.com/thread/223790-how-to-set-the-basemap-layers-of-a-webmap-using-the-arcgis-python-api


texts: CORS Printing Issue with ArcGIS Javascript

texts:
I'm developing an application using the ArcGIS Javascript API and everything has been working fine until I've tried to incorporate the print widget using printtask.  When only the base layers are accessed in the code, the map is printed just fine.  But when I include any layer running on my local instance of ArcGIS Server 10.2 (dynamic or feature layers) the print operation fails.  I assume this is because I need a "proxy" file although I thought I read it's possible to utilize CORS.  I'm using the latest FireFox browser, ArcGIS Server 10.2 and I have configured my code to include:
esri.config.defaults.io.corsEnabledServers.push("10.111.123.234:6080");
esri.config.defaults.io.corsEnabledServers.push("sampleserver1a.arcgisonline.com");
esri.config.defaults.io.corsEnabledServers.push("sampleserver6.arcgisonline.com");
esri.config.defaults.io.corsEnabledServers.push("tasks.arcgisonline.com");
esri.config.defaults.io.corsEnabledServers.push("server.arcgisonline.com");

The IP address above is my local PC pointing to the ArcGIS Server port of my services.  I have tried just localhost:6080 as well.

Either way, the printing functionality still fails so I'm not sure if CORS is valid for this operation and that I have to use a proxy file anyway?


texts: Adding a field with a unique value to vector layer in kml file (need this to burn values with gdal.RasterizeLayer)

texts:
I would like to create a field ID in my kml file, with a unique value so that I can use gdal.RasterizeLayer to burn the value to an image like this:

gdal.RasterizeLayer(target_ds, [1], source_layer, options = ["ATTRIBUTE=ID"])

This means I need to create an attribute called ID in my .kml file so that gdal.RasterizeLayer will be able to recognize this and burn the correct ID value. I tried this and it seems like I create a new field but when I run gdal.RasterizeLayer, nothing is rasterized. My current kml file has several Placemark/Polygon objects and I wish to give each of them a unique id.

import pykml.parser as kml_parser
from lxml import etree

kml_file = 'test.kml'

with open(kml_file) as f:
    root = kml_file.parse(f).getroot()

for i, Placemark in enumerate(root.Document.Folder.Placemark):
    Placemark.addattr('id', i + 1)

tree = etree.ElementTree(root)
tree.write('new.kml')

Does anyone know what the proper way to create attributes in kml file is?


texts: Altering font size using label expression?

texts:
I have two label classes displaying dynamically over my parcel polygons. The owner info is the major label with larger font size and bellow it I am displaying some property type info in smaller font.

The problem is that the 2 label classes just will not display the way I want them to; Major label at the top, minor label just bellow it.  I've used maplex and the regular labeling engine and played with the options but it just does not do what I want and the labels jump all over the parcels.

The only way that works is to export the labels into a feature linked annotation size then adjust the label positions manually.  This is OK but i need to change the scale often on this map.

Can I write the label expression so that I can have 2 different fonts and sizes within one label?  (Something similar to using text tags with text graphics)