texts: Use R to combine information in a raster and a point vector according to a separate polygon vector associated with each

texts:
This builds off a previous question I had about counting the pixels in a raster that contains at least 1 point (Use R to count the number of pixels in a raster file for which there are at least one point object present) which was perfectly answered by @dbaston.
I guess my real problem is a bit more complex and I am curious how to arrive at a similar solution with an added layer of complexity, where instead of just counting up the pixels with points for the entire raster, the counts of pixels with points is made separately for a number of uniquely ID'ed polygons.
Here is a visual representation of the problem

Here is the ideal outcome

Reproducible Example
library(raster)
set.seed(123)
#===========================
# POINTS
#===========================
# create points
points <- structure(list(longitude = runif(10), latitude =runif(10)), .Names = c("longitude", "latitude"), class = "data.frame", row.names = c(NA, -10L))

# set coordinates for spatial data frame
coordinates(points) <- cbind(points$longitude , points$latitude)

# assign a CRS to spatial object
proj4string(points) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#===========================
# RASTER WITH DATA
#===========================
# create raster
rast <- raster(xmn=0, xmx=1, ymn=0, ymx=1, res=0.25)

# add data
rast$data <- rnorm(length(rast), mean = 10, sd=1)

#===========================
# POLYGONS
#===========================
square <- rbind(c(1, 0.5, 0.5, 0.5, 0.5, 1,
                  1, 1, 1, 0.5, 1, 0.5),
                c(0.5, 0, 0, 0, 0, 0.5,
                  0.5, 0.5, 0.5, 0, 0.5, 0))

ID <- c("Poly1", "Poly2")

polys <- SpatialPolygons(list(
  Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=TRUE))), ID[1]),
  Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=TRUE))), ID[2])
))

#===========================
# VISUALIZE
#===========================
plot(rast)
plot(points, pch=1, add=TRUE)
plot(polys, lwd=3, add=TRUE)

That is where I am stuck. I think this should be easy but it is taking a fair bit of time to figure out for someone who is trying to learn how to work with spatial data in R. Thanks in advance for any insight you can share!


texts: Upload GML to Geoserver as Layer through REST

texts:
I have to upload a GML file to my geoserver through REST. I tried to use importer but it doesn't work (even manually). More specifically, the importer module never generate tasks automatically for GML. When I use another REST call to generate tasks, it returns max-stack-limit error. Further, I noticed that GML is not in the supporting list of importer extension.

Here is the code I used to create a import event:

          var uri = 'http://192.168.0.30:8080/geoserver/rest/imports';
          var data = JSON.stringify({"import": {
                 "targetWorkspace": {
                    "workspace": {
                       "name": "wfsresult"
                    }
                 },
                 "targetStore": {
                    "dataStore": {
                       "name": "WFSSourceResult"
                    }
                 },
                 "data": {
                   "type": "file",
                   "file": "/home/zhangj/geoserverDownload/data/wfsresult/u4f3809a0-56a3-11e7-894d-8dc4d5efd317/result.xml.kml"
                 }
              }});
          Request({
            uri: uri,
            method: 'POST',
            auth: auth,
            headers: {
              'content-type': 'application/json'
            },
            body: data
          }, (err, response, body) => {
            if(err) {
              reply(err);
            } else {
              reply(body);
            }
          });
        }

Am I using importer extension correctly?

If importer doesn't fit my situation, what should I use?

I got that wfs transaction may be a good solution. How should I use it?


texts: Selecting a series of grouped layers within a specific extent

texts:
Is there a way to select grouped layers within a specific extent with ArcGIS 10.5?

For example I have a constraint mapping, with various thematic layers grouped by interest. I would like to know which one are applicable in my project area, without having to select and clip each layer individually (too many of them). Ideally a tool that would let me grab the group of layer and create a shp or layer of features present within my work area would be great.


texts: Set multiple values to nodata - Erdas Imagine 2015

texts:
I am having trouble with the data ranges in the Landsat 8 surface reflectance product which I downloaded from USGS. The Lansat 8 surface reflectance handbook says the valid range is 0 - 10,000 range but the product may come with values between -2000 to 16000. I am using the level 2 landsat 8 surface reflectance product with a view to detecting sediment plumes exiting glaciers by thresholding.

I want to know how to do the following on Erdas Imagine 2015 to:

exclude the values outwith the valid range (0 - 10,000)
apply the scale factor of 0.0001 in order to convert reflectance values to a scale from 0 to 1.

Is there a way to apply the Set to NoData function in the Spatial Model Maker to multiple values - say those values that are < 0 and > 10,000 ?


texts: "Invalid Data Source" When trying to open many w001001.adf files

texts:
I have a folder with more than 9,000 subfolders, each with ESRI GRID files. Here's the setup of the directories and subdirectories.

Data-folder
- be_rasters
-- 692500_725000
--- dblbnd.adf
--- hdr.adf
--- metadata.xml
--- prj.adf
--- sta.adf
--- w001001.adf
--- w001001x.adf
-- (thousands more folders with this naming convention format and files within them named the same thing)

I'm trying to work in QGIS 3.8 to open the 9000+ w001001.adf files in the folders in be_rasters directory. In QGIS, I go to Layers >> Add Vector Layer, then choose "Directory" as a Source Type and "Arc/Binary Coverage" in the "Source" section. Then under "Vector Dataset(s)" I select be_rasters.

But when I click "Add," I get this error.

Invalid Data Source: path/to/Data-folder/be_rasters  is not a valid or recognized data source.

How do I load all the w001001.adf files into QGIS?


texts: PostgreSQL : Question about dropping tables inside funtions

texts:
I am writing some PostgreSQL code where I have a function which is used multiple times in a single query with different input values. Inside the function I create some tables and use them.

What is the correct way to handle those tables? should I drop them at the end of the function (right before the RETURN), or should I use temp tables?


texts: Overlaying a .png image that comes with a world file in UTM with OpenStreetMap data?

texts:
I have received a geocoded .png image (2380 x 2396 pixels) that comes with a world file containing georeferencing information in UTM format:

worldFileParams = [0.5,0.0,0.0,-0.5,690920.0,5335407.5]

Now, I would like to overlay this image with map data. I thought about using OSM data for that purpose. In order to solve the task, I tried the following workflow:

Download shapefile of OpenStreetMap data from https://mapzen.com/data/metro-extracts/
Load both shapefile(s) and .png raster image in QGIS.

However, this surely naive approach hasn't worked, as both datasets are not co-located in QGIS, which is probably caused by the fact that QGIS data is provided in WGS84 (EPSG:4326) coordinates, while the png world file data is in UTM32N (EPSG:32632). When I try to reproject the png file, I receive the error

  ERROR 1: latitude or longitude exceeded limits
  ERROR 1: Too many points (2401 out of 2401) failed to transform,
  unable to compute output bounds.

Thus, long story short: What can I do to reach my goal? Although using GDAL, QGIS, MATLAB and OpenStreetMap would come in handy for me, I am not bound to using these things.


texts:
I'm trying to set my shields on the map which contain only image (I don't want any text at all - only img). I can't change size of img.
I tried markers too but I can't get them orientated like shields (no rotation, shield on the line).

How can I achieve that?

As you can see icons are way too big but I don't want to change size of files.

    #rowerowe {
    [colour='red']{
    shield-file:  url("D:\Leaflet\img\red.png");

    shield-face-name:"Arial Bold";
    shield-size:16;
    shield-fill: black;
    shield-line-spacing:5;
    shield-min-distance:1;
    shield-allow-overlap:false;
    shield-placement:line;
    shield-dx:150;
    shield-spacing:450;
    }
 }

Here's the same with markers, I can change size but can't change orientation to "top".

#rowerowe {
    [colour='red']{
    marker-file:  url("D:\Leaflet\img\red.png");
    marker-width:24;
    marker-height:24;
    marker-placement:line;
    marker-spacing: 200;

    }
 }


texts: Remove reference to raster layer object PyQGIS

texts:
I am developing a plugin to load a set of raster files on QGIS. When there is a large number of files to be loaded, the main UI freezes. So I use QgsTask to move the creation of raster layers to a thread. This thread then returns a list of raster layer objects. I was able to load these layers on QGIS.
When I remove the layers from canvas by using QGIS UI or QgsProject removeMapLayers, the layers are removed. However, it seems that the raster layer objects still exist. I cannot move or remove the source raster files on disk. It says that the files are still being used by qgis.bin. This did not happen when I did not use the thread.
Here is my code:
Create raster layers task:
task = QgsTask.fromFunction('Load images', createRasterLayers, params = [result["files"]], on_finished=self.load_layers)
task.progressChanged.connect(lambda: self.dlg.he_progress.setValue(task.progress()))
task.taskTerminated.connect(lambda: self.dlg.he_progress.reset())
QgsApplication.taskManager().addTask(task)

The function that creates a list of raster layers, but not load yet:
def createRasterLayers(task, params):

    # extract parameters
    images = params[0]

    n_photos = len(images)

    # set task progress
    task.setProgress(1)

    # then reload the photos
    count = 0
    n_processed = 0
    rlayers = list()
    rt = QgsRasterTransparency()
    rt.initializeTransparentPixelList(0,0,0)
    for img in images:
        fileInfo = QFileInfo(img)
        path = fileInfo.filePath()
        basename = fileInfo.baseName()
        r = QgsRasterLayer(path, basename)
        crs = r.crs()
        crs.createFromId(4326)
        r.setCrs(crs)
        # set transparent for the image
        r.renderer().setRasterTransparency(rt)

        if r.isValid():
            rlayers.append(r)
            count = count + 1

        # set progress
        n_processed = n_processed + 1
        percent = float(n_processed/n_photos) * 100
        task.setProgress(percent)

        if task.isCanceled():
            raise Exception('Task canceled!')

    return {'rlayers': rlayers, 'task': task.description()}

And finally, the function that loads the layers on map:
def load_layers(self, exception, result=None):

    if exception is None:
        for r in result['rlayers']:
            QgsProject.instance().addMapLayer(r, addToLegend=True)


texts: How to get status of a layer in geoserver using gsconfig?

texts:
I am getting all published layers from geoserver using gsconfig using following code:

from geoserver.catalog import Catalog
cat=Catalog("http:// localhost:8080/geoserver/rest")
myResource = cat.get_resources(store='storeName', workspace='workspaceName')
layers = cat.get_layers(myResource).

Now I want to get all layers which are not published stored (yet to be published) or status of a layer which means published or yet to be published . How can I?


texts: Spatial weight matrix with boundary effects in R?

texts:
I need some help to finalize code for a spatial weight matrix that uses a nearest neighbor definition within regions, but does not allow neighbors from across a political border. For example, physically adjacent counties are neighbors, but physically adjacent counties separated by a state boundary are not neighbors.

I am working with an ESRI shapefile, with about 40,000 polygons.

The nearest neighbor bit is straightforward using functions from spdep. And I have collected a code example that creates a boundary effect:

# xxx is the name of the shape file

unique(xxx$country)
length(unique(xxx$countrycode))

# B_domestic = Indicator function taking value 1 if two regions belong to  same          country and 0 if not

B_domestic <- bdiag(
        matrix(1,length(xxx$countrycode[xxx$countrycode=="001"]),
            length(xxx$countrycode[xxx$countrycode=="001"]))
        matrix(1,length(xxx$countrycode[xxx$countrycode=="002"]),
            length(xxx$countrycode[xxx$countrycode=="002"])),

##  and so on .  . .
##  then check the object generated

dim(B_domestic)
class(B_domestic)

##  convert it to a matrix
B_domestic <- as.matrix(B_domestic)

## then the matrix is converted to a list W object

B_domestic_listW <- mat2listw(B_domestic); class(B_domestic_listW)

## and back to a weights object (I think? documentation is pretty sparse)

B_domestic_2 <- as_dgRMatrix_listw(B_domestic_listW)

It seems that the something from this process then becomes an argument for creating the nearest neighbors, but what is the cleanest way to do this? And in which function?

Using spdep, I have:

data4NN <-  knn2nb(knearneigh(coords, k=4), row.names = ID)
##row-standardized

W.list <- nb2listw(data4NN,style="W")
summary(W.list)

where coords are the previously generated centroid coordinates of my original polygon file, and ID is a unique identifier in the polygon file; 4 nearest polygons are neighbors.

Object W.list then becomes the object used in my spatial regression specification (SAR,SD,SEM)

I feel like I have the pieces, but don't understand how they relate to one another to get the outcome I'm looking for. Specifically, where does the argument to limit the adjacency by the political boundary fit into generating a nearest-neighbors spatial weight matrix?


texts: Calling a AddIn from another class or AddIn?

texts:
I've tried several ways but I can't find anything working. Is there any possibility to call an AddIn from outside? I want to fill a AddIn Combobox when the OnClick Event of an AddIn Button is fired, but don't know how to address the Combobox.

Have you any Ideas?


texts: Add attributes to an output shapefile schema (via meta)

texts:
I am trying to copy the schema of an existing shapefile and add to it in an output shapefile. Somehow I am not getting anywhere with this...though the docs seem very clear (and the download page : https://pypi.python.org/pypi/Fiona) and it is just adding to a dictionary.

For this object, I have this schema

      c = fiona.open(r'I:\It_24\115507_Road_Shields_Label_Processing\data\test_data.shp')

      c.schema
      {'geometry': 'LineString', 'properties': OrderedDict([(u'ID', 'float:11'), (u'OSM_ID', 'float:19'), (u'NAME', 'str:254'), (u'TYPE', 'str:254'), (u'TUNNEL', 'int:6'), (u'BRIDGE', 'int:6'), (u'ONEWAY', 'int:6'), (u'REF', 'str:254'), (u'Z_ORDER', 'int:6')])}

and I want to return:

      c = fiona.open(r'I:\It_24\115507_Road_Shields_Label_Processing\data\test_data_copy.shp')

      c.schema
      {'geometry': 'LineString', 'properties': OrderedDict([(u'ID', 'float:11'), (u'OSM_ID', 'float:19'), (u'NAME', 'str:254'), (u'TYPE', 'str:254'), (u'TUNNEL', 'int:6'), (u'BRIDGE', 'int:6'), (u'ONEWAY', 'int:6'), (u'REF', 'str:254'), (u'Z_ORDER', 'int:6'), (u'shield_type', 'str:254'), (u'label', 'str:254'),  (u'label_len', 'int:10'), (u'zoom', 'int:10')])}

to do so I am running:

   def process_file(self, inFile, outFile):
        with fiona.open(inFile, 'r') as input:
            meta = input.meta
            # create new fields for the new schema
            meta['schema']['properties']['shield_type'.encode("utf-8")] = 'str:254'
            meta['schema']['properties']['label'.encode("utf-8")] = 'str:254'
            meta['schema']['properties']['label_len'] = 'int:10'
            meta['schema']['properties']['zoom'] = 'int:10'

            with fiona.open(outFile, 'w', **meta) as output:
                  for item in input:
                    n = item.copy()
                    new_data_attributes = function_using_some_regex_parsing(item['properties']['REF']
                    ...

Right now running this I just get an empty dictionary...I'm not sure what I'm missing...

      c = fiona.open(r'I:\It_24\115507_Road_Shields_Label_Processing\data\test_data_copy.shp')

      c.schema

      {'geometry': 'LineString', 'properties': OrderedDict()}


texts: Can't change render type from single band color data in QGIS 3.16 on DSM LiDAR layer

texts:
I can't change the render type from single band color data in QGIS 3.16 on a DSM LiDAR layer.
When I click on the drop down arrow by render type there are no other options. I would like to add multiple colors to the DSM layer.


texts: Travel Time Isochrone Map - Http Post request

texts:
I am trying to send a request to Travel Time Platform for an isochrone using the API - https://www.traveltimeplatform.com/

// Coordinates for Manchester Central Foodbank in Lat Lng
    var ManchesterCentralFoodbank = "53.464774, -2.23179";

    // The departure time in an ISO format
    var departureTime = "2019-12-04T09:00:00-0500";
    // Travel time in seconds. We want 1 hour travel time so it is 60 minutes x 60 seconds.
    var travelTime = 60 * 60;

    const url = "https://api.traveltimeapp.com/v4/time-map";

    // API key and ID
    var APPLICATION_ID = "hidden_here";
    var API_KEY = "hidden_here";

    var request = {
        departure_searches: [{
        id: "starting_location",
        coords: ManchesterCentralFoodbank,
        transportation: {
            type: "public_transport", },

            departure_time: departureTime,
            travel_time: travelTime, }],

    arrival_searches: [] };

    // This function sends the request and draws the result onto the map
    function sendIsochrone(url, drawIsochrone) {

    let xhttp = new XMLHttpRequest();
        // event listener for when the HTTP state changes
        xhttp.onreadystatechange = function () {

        //a successful HTTP request returns a state of DONE and a status of 200
        if (xhttp.readyState === XMLHttpRequest.DONE && xhttp.status === 200) {
        drawIsochrone(JSON.parse(xhttp.responseText));
        }
      };

    // the request itself
    xhttp.open("POST", url)
    xhttp.setRequestHeader("X-Application-Id", APPLICATION_ID);
    xhttp.setRequestHeader("X-Api-Key", API_KEY);
    xhttp.setRequestHeader("Content-Type", "application/json; charset=UTF-8");
    xhttp.responseType = "json";

    // passes the information stated in the request variable
    xhttp.send(JSON.stringify(request));

    }

    function drawIsochrone(data) {

        console.log(data)

    }

I keep getting a CORS error message  -  "Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at file"
-   I cant see any problems with my code and I am not sure where I'm going wrong.

I tried making these changes to my firefox browser https://support.mozilla.org/en-US/questions/1264280 but to no avail.


texts: QGIS Joined Layer not displaying points visible on pre-join layer (No Extent)

texts:
I have a shp file which is the output of the Toolbox Check validity tool in QGIS 3.4.3
The points display for me correctly but if I then join the layer with a csv attribute using the toolbox Join Attributes by Field Value to create a new layer the points no longer display on my map.

If I then select 1 of the features using the attribute table and use the zoom map to selected feature I get the error

  Cannot zoom to selected feature No extent could be determined

I have also tried to associate the attribute table using the Layer properties Join but I end up in a similar position.

Looking at the layer properties it seems the extent value is lost in the processing but I'm not able to restore it using the "Update Extent" option on the Joined layer properties.

Anyone know how to solve my issue?

Updates

I tried the multipart to single part suggested in the comments but
that didn't work.
Also if I try to filter the valid layer itself manually by creating a valid expression to only show the points with a specific value ColName in ('ID1') (which I know in the csv file) all the points disappear and I cannot zoom to the selected. This seems to me to point to an issue with my points layer not the join but not sure what.

Joined Layer

Encoding windows-1258
Geometry Point (MultiPoint)
CRS EPSG:4326 - WGS   84 - Geographic
Extent Empty
Unit degrees
Feature count 46

Valid Layer

Encoding windows-1258
Geometry Point (MultiPoint)
CRS EPSG:4326 - WGS 84 - Geographic
Extent -10.4506078999999996,51.4353689000000003 -6.0394994000000004,55.3615765999999994
Unit degrees
Feature count 4,735


texts: Distance between two points on QGIS 3.4.7 and Python 3.7

texts:
I am under the 2154 CRS (RGF 93).
I want to calculate the distance in meters between two point knowing their coordinates.

Let's consider two points : A(474828.85, 6756169.31) and B(874895.75, 6756159.5)
If I measure the distance manualy I get 69.6m.

I want to get this result through the Python console. I have read two ways of doing it: QGIS Python Distance between points

Nevertheless, it fails on my case.

>>>point1 = QgsPoint(474828.85, 6756169.31)
>>>point2 = QgsPoint(874895.75, 6756159.5)
>>>distance = QgsDistanceArea()

This is ok, but the next line sends an error message:

>>>m = distance.measureLine(point1, point2)
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
TypeError: QgsDistanceArea.measureLine(): arguments did not match any overloaded call:
  overload 1: argument 1 has unexpected type 'QgsPoint'
  overload 2: argument 1 has unexpected type 'QgsPoint'

The other method presented in the link fails too:

My question is : "How can I get the distance given in meters ?"


texts: Project small raster without causing distortions in R

texts:
I want to change the units of small raster files from degrees to meters using R.
Is there any way to do this without projecting and hence altering the extent and shape of the rasters in the process? Or is there a projection that doesn't distort small areas at all?
I am aware that projections always distort the area of interest to some degree and that this is their desired property. However, my rasters are very small (less than 10x10 km) and have a very coarse resolution (ca. 300x300 pixels), and hence I would prefer to simply change my unit of measurement from degree to meters without altering the raster at all (and just accept that this causes some minor geographic inaccuracies)
I've tried both using the appropriate UTM zone or a Lambert Azimuthal Equal Area projection centered on my raster, but both projections slightly change the extent of the raster.


texts: Compare slopes between adjoining line segments

texts: Using a DEM, I assigned slope values to individual segments in a street network feature class. In some areas the slope values seem suspect, such as where a highway overpasses a local road. Is there a method to find those segments whose slope is radically different from adjoining segments (e.g. 5 to 25%) to locate possible inaccuaracies?


texts: How to install Proj4js.js in my localhost directory?

texts:
I'm developing a GIS website for vineyard management. I am using Openlayers with the WMTS IGN (French institute for geography) datas. It works fine.

I also would like to add the possibility for the user to switch to a geology map provided by the BRGM (geological french institute) geoservice.

The issue is that IGN datas avec in Webmercator projection (EPSG:3857) and BRGM datas avec in EPSG:27572 (NTF (Paris)) and i would like to keep map center and bounds while switching from one to the other one. After reading a lot on the internet, i found i needed to use Proj4js for projection transformation. I also have vestor shapes (vineyard fields) and i would need to visualize in both projections by transforming them.

I am quite used to implement js librairies (by copying js files in my localhost and include it in the header of my php/html file) but in the case of proj4js, it looks different. The gituhub page only gives this information for installation :

npm install proj4
bower install proj4
jam install proj4
component install proj4js/proj4js

I have no idea how it works and i can't get any clear information from what i read. In the github files, there is no proj4js.js file... Is there somewhere a clear tutorial for installation ? Does someone have experience on this installation ?

I tried a "classical" installation downloading 1.1.0 from this link http://download.osgeo.org/proj4js/proj4js-1.1.0.zip but it seems to be an old version and i'm getting errors with the test code below :
    "lcc:forward: llInputOutOfRange: NaN : NaN"

Here is my test code :

<head>

<script src="js/proj4js/lib/proj4js-compressed.js"></script>
<script src="js/proj4js/lib/defs/EPSG3857.js"></script>
<script src="js/proj4js/lib/defs/EPSG27572.js"></script>
<script src="http://openlayers.org/api/2.13.1/OpenLayers.js"></script>

</head>

<body>

<script type="text/javascript">

var source = new Proj4js.Proj("EPSG:3857");
var dest = new Proj4js.Proj("EPSG:27572");

Proj4js.reportError = function (msg) {
    throw new Error('Whoops: ' + msg);
};

var position1 = new OpenLayers.LonLat(0, 5700000);

    console.log(source);
    console.log(dest);
    console.log(position1);

    console.log(Proj4js.transform(source, dest, position1));

</script>

</body>


texts: Does ArcGIS Desktop automatically install Python when there is update?

texts:
We noticed an office-wide installation of Python 2.7.15 this morning at 5:30am.

The installation pushed Python 2.7.15 into a standalone folder on the C:/ drive of all our PCs.  If there are two instances of Python on the machine, ArcGIS Desktop cannot recognize which Python instance to use, thus dropping functionality with the ArcPy module.

This is easily tested by launching the arcpy window in ArcMap and typing in import arcpy.  The application crashes...

I assumed it was an auto-install setting in ArcGIS, but then we would expect the install to occur within the ArcGIS installation folder/directory and not a standalone install outside of programfiles.

Has anyone else experienced this?

We are on a campus that has fairly strict firewalls and the IT department was unaware of the automatic install.

This seems like a bug that could be reported but wanted to reach out to this community first.


texts: ODBC connection to ArcSDE table shows old version in ArcMap?

texts:
I use an ODBC connection to an ArcSDE table to return a current list of features for reporting purposes.  The connection is pointed at a Microsoft SQL ArcSDE server and the connection info matches what I see in ArcCatalog.

The problem is that the ODBC connection returns a version of the table that is several months old - it shows some features that have been removed and doesn't show some new features.

Using the same connection information in ArcMap, I see a current version of the table and the current set of features.

I have created a new ODBC connection, tried refreshing the connection, and tried using the connection in Microsoft Access and R.

Any suggestions?

This question referred to a versioned table.  As it turns out, an editor had checked out a parent version of the ArcSDE some time ago and never checked it back in.  This was preventing subsequent child versions from being compressed back to the DEFAULT table. Child edits were being kept in a delta table but still showing in ArcMap; the ODBC connection wasn't aware of the delta table and was showing DEFAULT as it existed when the rogue parent version was checked out.


texts: Find polygon that contains all linestring records in PostGIS table

texts: I'm still pretty new to this stuff, so bear with me.  In PostGIS how could I find a polygon that would contain all the linestring records in a given table?


texts: Geoserver Explorer cannot connect to Geoserver returning error 404

texts:
This is to share my experience on a problem I had and how I solved it:

I have just installed GeoServer and it works fine, however, I wanted to give a try to the GeoServer Explorer but it could not connect to GeoServer. I am not sure what I have done wrong.

I am using Windows 10 and QGIS 2.18.4 . I get the following error when trying to create a new catalog from QGIS:

Tried to make a GET request to http://localhost:8080/Geoserver/Web/rest/workspaces.xml but got a 404 status code:
HTTP ERROR: 404
Problem accessing /Geoserver/Web/rest/workspaces.xml. Reason:
    Not Found

Powered by Jetty://
Traceback (most recent call last):
  File "C:/Users/mygis/.qgis2/python/plugins\geoserverexplorer\gui\explorer.py", line 101, in run
    command(*params)
  File "C:/Users/mygis/.qgis2/python/plugins\geoserverexplorer\gui\gsexploreritems.py", line 551, in populate
    raise e
FailedRequestError: Tried to make a GET request to http://localhost:8080/Geoserver/Web/rest/workspaces.xml but got a 404 status code:

HTTP ERROR: 404
Problem accessing /Geoserver/Web/rest/workspaces.xml. Reason:

Not Found


texts: PostGIS - raster coodinates - no result for negative values of GPS

texts:
I have data in raster and try to run this query

WITH ref As ( SELECT (ST_SetSRID(ST_Point(-20.12, 51.50), 4326)) AS pt )
SELECT ST_Value(raw_data, 1, ref.pt) AS val_0
FROM model_data

The problem is, that for negative longitude (-20.12) I do not get any output (latitude can be +- 51.50 and I still getting results).

My data st_upperleftx = -0.125 and st_upperlefty = 90.125. I can add offsets to coordinates lat/lon manually, but that is not what I want, because every raster can have a different upperleft corner. I thought, that postgis should take care of this automatically, but I am obviously missing something somewhere.


texts: Select highest value in multiple polygons

texts:
I have a polygon shapefile and a point shapefile. The point shapefile has the following attributes: X, Y,height (derivated from a reaster DEM) and the polygon unique ID.

I want to select the highest point on each polygon, but I don't know how to use the correct expression to accomplish it. Can anyone help me?


texts: QGIS crashing when trying to open PostGIS layer

texts:
I'm trying to open a fairly large dataset from a postgres database in QGIS. The table contains over 100 million rows, and it about 60gb in size.

When I try to open this table in QGIS (have tried 2.8 and 2.14) it instantly crashes. It doesn't even try to open it. QGIS just locks up completely, no memory usage, no CPU or network activity.

I've tried zooming in to a small area before trying to load the data, but it does the same thing. Also, similar, smaller datasets open fine in the same way.

I was wondering if anyone has had a similar issue? Do we know if there is a limit to the size of a dataset that QGIS will handle from postgres?


texts: QGIS and LRS. Define points along line from distance in .txt file

texts:
I'm working with QGIS 2.14.3. and am a fairly novice user.

I have a shapefile layer with polylines, each line is identified by a unique 'LCS_CODE' field.
In a spreadsheet (imported as a text file layer) I have speed sign data which I want to display. There are three columns ('SPEED','LCS_CODE','DISTANCE'), the 'SPEED' of the sign, the 'LCS_CODE' it is on and the 'DISTANCE' from the start of the LCS_CODE it is located. There may be more than 1 speed sign per LCS_CODE. Other than the distance field there is no geographic info.
I want to get QGIS to display the speed sign at the specified distance along the LCS_CODE line.

My first thought was to use LRS but I cannot get it to 'see' the text file layer to load it as the point layer.
I then created a points layer with the same 3 fields (SPEED, LCS_CODE, DISTANCE) with the thought that I would copy/paste the spreadsheet data straight into the attribute table but cannot see a way to do this without artificially creating points for each row first.

How else might I achieve this, to me it sounds reasonably simple however it's obviously beyond me!!


texts: UpdateCursor for field calculation based on unique values

texts:
UPDATE TO PREVIOUS QUESTION: I have an address feature class. I want to use an UpdateCursor that focuses on two fields: StreetName and Block_ID. For every instance of a unique StreetName, I want the Block_ID to start at 10 and add multiples of 10. For example:

StreetName    Block_ID
Main St       10
Main St       20
Main St       30
1st St        10
1st St        20
3rd St        10
3rd St        20

I have the following code but I know it is way off.....I was just testing that it was populating the Block_ID:

def unique_values(table , field):
    with arcpy.da.SearchCursor(table, [field]) as cursor:
    return sorted({row[0] for row in cursor})
# myValues is the list of unique street names identified in SearchCursor
myValues = unique_values(Table, 'STREET')
arcpy.AddMessage(myValues)

# Here is where we calc block id
strtFld = 'STREET'
blkFld = 'BLK_NO'
update = Table
i=0

with  arcpy.da.UpdateCursor(update, [blkFld, strtFld]) as cursor:
    for row in cursor:
        row[1] = myValues[i]
        row.setValue(row[0], i=i*10)
        cursor.updateRow(row)


texts: How to delete vertices with the same direction

texts:
I have feature with a big amount of vertices. But these vertices are totally useless and make Qgis editing too slow. Lets look on image with example:

These vertices inside of black line I want to delete. But I want to preserve shape of red line as well. So I need to let vertex outside of black lines.


texts: How to get the file name of current style using PyQGIS?

texts:
I get the name of the current style file using layer.styleURI().  Then, I change the style file using layer.loadNamedStyle(different_filename).

However, when I call layer.styleURI() again, it returns the same name as in previous call.  How can I get the file name of the current style in QGIS using Python?


texts:
In pyQGIS i tried this code (contains also debug prints):

def joinLayers(inputlayer, layer, index_field, name = "JOINED", dataProvider = "ESRI Shapefile"):
 a = []
 b = []
 outpath = getPath(inputlayer.source(), name)
 inputfields = inputlayer.pendingFields()
 nameoflayer = layer.name()
 toappend = layer.pendingFields()
 for fild in toappend:
    a.append(fild.name() + "j")
    fname = fild.name()
    newfield = QgsField(unicode(fname+ "j"), fild.type(), fild.typeName(), fild.length(), fild.precision(), fild.comment())
    b.append(newfield.name())
    inputfields.append(newfield, QgsFields.OriginJoin)
 print a
 print b
 output = QgsVectorFileWriter(outpath, dataProvider, inputfields, inputlayer.crs(), layer.wkbType(), inputlayer.dataProvider().encoding())
 del output

Main problem is renaming fields in for loop over QgsFields.

a

[u'UWIj', u'Wellj', u'UBHIj', u'Boreholej', u'Statej', u'Jurisdictij', u'Offshorej', u'Basinj', u'Sub-basinj', u'Latitudej', u'Longitudej', u'Legislatioj', u'Titlej', u'Operatorj', u'Classj', u'Structurej', u'Seismic loj', u'Drilled byj', u'Rigsj', u'Spud/Kickj', u'TD Datej', u'Rig releasj', u'Start deptj', u"Driller'sj", u'Oil showj', u'Gas showj', u'Condensatej', u'Oil fieldj', u'Discoveryj', u'Statusj', u'Well headj', u'Basic relej', u'WCRj']

b

[u'UWIj', u'Wellj', u'UBHIj', u'Boreholej', u'Statej', u'Jurisdictij', u'Offshorej', u'Basinj', u'Sub-basinj', u'Latitudej', u'Longitudej', u'Legislatioj', u'Titlej', u'Operatorj', u'Classj', u'Structurej', u'Seismic loj', u'Drilled byj', u'Rigsj', u'Spud/Kickj', u'TD Datej', u'Rig releasj', u'Start deptj', u"Driller'sj", u'Oil showj', u'Gas showj', u'Condensatej', u'Oil fieldj', u'Discoveryj', u'Statusj', u'Well headj', u'Basic relej', u'WCRj']

After write this fields into QgsVectorFileWriter, some fields doesn´t have "j" suffix.

For example attribute Condensate, see attribute table in qgis:


texts: Aggregating census level data to county level in QGIS?

texts:
I'm trying to take census tract linked data (in the form of a CSV I joined to a census tract level map of Washington State) and create a county-level map with it. My census data also has the county listed as a separate variable for each census tract (so basically what county each census tract belongs to). How can I take these census data and aggregate them to a county level in QGIS? More specifically, I'm trying to aggregate median income and percent of unemployment from the census level to the county level.

I've tried spatially joining the census map with the data to a county level map however the new values in the attribute table do not appear to be correct.


texts:
I've received a .dxf file of a small area in Sicily. Now I have a problem with the coordinate system. It's not the usual European UTM or similar. I absolutely cannot find what coordinate system it is.
The extent is
-4011.0709999999999127,-4410.8940000000002328
-790.0000000000000000,-2388.5529999999998836
I have been using QGIS (without Python) for a long time and I know my way around. Here I have not found a way to get the EPSG from the metadata. or is there a way?
I googled a lot but couldn't find the system.


texts: Python Geopandas reading shapefile taking too long

texts:
I'm relatively new to Python and GIS.

I am attempting to read in a 1.1 GB shapefile using geopandas. I let the code run for an hour and it is still not complete. Is this normal? How can I optimize the time it takes for the file to be read?

Here's my code:

polygons = gpd.read_file(r"path_to_shapefile.shp")

Edit: My system is an SSD, 64bit Python, and 8gb of RAM.

Edit 2: Not sure what was different, but running this code a third time, the shapefile was able to be read by Geopandas in less than 5 minutes.


texts: How to send variables from gis addin own web browser to gis file

texts:
I am a complete beginner with GIS.

So I want to be able to go to the coordinates in arcmap I fill in in my php page. I've got this now:

The variables inserted in the textfield should be the coordinates where I want to go to. The text inserted in X: and Y: are called $x and $y. How do I accomplish this?


texts:
I have a simple OpenLayers map, where I would like to display the OSM map (the base map, in EPSG:900913) with some points on top of it (a vector layer, in EPSG:3763).

I use proj4js. I do convert a point from EPSG:3763 to EPSG:900913, to set the map center. I also added the MousePosition control, and I can see the coordinates in the displayProjection system. so, the conversion runs fine between the two EPSGs.

I created a new vector layer, with the EPSG:3763 projection. I add a point to that layer pragmatically, in the layer projection EPSG:3763.

The point in the layer goes to Africa, into the ocean, around 0,0.

My questions is:
Does OpenLayers can reproject the entire layer? If it does, what am I doing wrong?

From the documentation, it not not clear if OL does the reprojection or not of the layer. In one example, the reprojection is done by MapServer. On the other, the reprojection is done when reading the geojson.

An alternative solution, is to convert the point coordinates to EPSG:900913 before adding them to the vector layer. But I really would like to see how to reproject an entire vector layer, in a different coordinate system.

My code is:

        <title>Multiple projections</title>
        <meta http-equiv="Content-Type" content="text/html; charset=UTF-8" />

        <link rel="stylesheet" type="text/css" href="http://dev.openlayers.org/releases/OpenLayers-2.13.1/theme/default/style.css" />

Happy mapping!


texts: For Loop Batch Copy Raster TIFF to BMP in ArcMap

texts:
I need to convert TIFFs to BMPs using CopyRaster from ascii to raster output TIFFs. I am not sure how to define the "outRaster" so that it will convert the rasters to BMP. I'm guessing that the line of code which I denote below needs additional consideration.

import os,arcpy,glob

rasterpath = r"C:\VMshared\small_example_valley2\snowdepthout"
outFolder = r"C:\VMshared\small_example_valley2\snowrast"
rasterList = glob.glob(rasterpath + "/*.tif")

for rasterFile in rasterList:
    outRaster = os.path.join(outFolder, rasterFile + "bmp") #THIS LINE??
    arcpy.CopyRaster_management(rasterFile,outRaster,"#","-3.402823e    +038","NONE","NONE","8_BIT_UNSIGNED","ScalePixelValue","NONE")

***NOTE - The code below works to convert to jpg or bmp although I would guess some syntax is missing.

import os,arcpy,glob

arcpy.env.workspace = r"C:\VMshared\small_example_valley2\snowtest2"
outFolder = r"C:\VMshared\small_example_valley2\snowtest3"

for rasterFile in arcpy.ListRasters("*.tif"):

    oName, oExt = os.path.splitext(rasterFile + ".jpg")
    outRaster = os.path.join(outFolder, oName + ".jpg")
    arcpy.CopyRaster_management(rasterFile,outRaster,"#","#","-9999","NONE","NONE","8_BIT_UNSIGNED","NONE","NONE")


texts:
I am drawing the habitat map within 500 Meters boundary. Once the works area changes, I need to update all the habitats (polygon) to align with the boundary (another polygon).

Recently, I am using editing shape tool to do it.

Is there any smarter way that I can align all polygons to the boundary at once?

I am using ArcMap 10.5.1


texts: looking for an esri compatible javascript engine

texts:
I work on an esri heavy team, i'm looking for a mapping engine that will give me flexibility to build on top of and leverage only the pieces I need and or add more as the app grows.  I've been using ESRI's API to great success for my indiviual code but leveraging their codebase into a reusable solution grows too large.

We have arcgis server 10.1 and create featureserver and mapserver services.  Along with those I will consume data thats created from various tabular locations and sent to my front-end via JSON/Ajax.

I've started on leaflet, and have began testing on esri-leaflet to some success, dynamic layers are drawing slowly, but we are currently upgrading to 10.1 so that may wash out.

Goal: find a lightweight javascript based mapping engine that will handle ArcGIS feature&mapserver services. Along with the ArcGIS services the engine needs to handle JSON data.  Geocoding with the service of my choice would also be great(bing, google, esri, or internally created)
Thx


texts: Add CSV file, join to feature class, field calculate and save output using ArcPy?

texts:
I am having some difficulty conceptualizing how the code should work for a custom ArcPy Script tool. Currently our work flow involves QA/QC and taking a csv file and joining it to a feature class. From there we calculate each field individually based on the values in the csv. This takes a tremendous amount of time. I am trying to automate this with a script. I am looking to create a tool allows one input (the csv file) and the tool will read the csv, join it to the feature class and then populates the fields automatically. I would also like to have some errors thrown if there is not able to join/calculate etc. I am not sure how this should flow.

import arcpy
import csv

def bringincsv(csv):
  with open(arcpy.GetParameterAsText(0), "rb") as csv:
   reader = csv.reader(csv, delimiter=",")
  next(reader):
  for row in reader:
   facilityid = str(row[0])
   diameter = double(row[1])
   hyperlink_cctv = string(row[2])
   hyperlink_rpt = string(row[3])

try:
 arcpy.env.workspace = 'CURRENT'
 inspection = arcpy.GetParameterAsText(0)
 arcpy.MakeTableView_management(in_table=inspection, out_view= 'Inspection')
 arcpy.AddJoin_management('Memphis.GIS.ssGravityMain', 'FACILITYID',
 'Inspection', facilityid, 'KEEP_COMMON')

arcpy.CalculateField_management('Memphis.GIS.ssGravityMain', diameter, ##The
csv files diameter)
arcpy.CalculateField_management('Memphis.GIS.ssGravityMain', hyperlink_cctv,
 ##The csv files values)
 arcpy.CalculateField_management('Memphis.GIS.ssGravityMain', hyperlink_rpt,
 ##The csv files values)


texts: How to export .mbtiles file from PostGIS to serve with S3

texts: We store our data in a PostGIS database. The database is updated roughly every week. However we do not want to serve vector tiles directly from our PostGIS but would much rather export a .mbtiles file to S3. The only option we see at the moment is exporting the PostGIS data to GeoJSON and afterwards using tippecanoe to bake the tiles. Is there a direct way to generate an .mbtiles file from a PostGIS database which can then be used by a static tile hosting service like S3?


texts: IO Performance of ArcPy in virtual machine

texts:
The following script just copies a file into RAM, performs some simple operations and writes the result back to hard disc:

import arcpy
import random
import time

arcpy.management.Delete("in_memory")
input_table = "C:\\data\\entrypoints\\table.shp"
memory_table = "in_memory\\table1"

# input from disc
t0 = time.clock()
arcpy.management.CopyRows(input_table, memory_table)
print time.clock() - t0

# some examplary operations and conversions
arcpy.management.AddField(memory_table, "NewField", "FLOAT")
tab = arcpy.da.TableToNumPyArray(memory_table, "*")
for row in tab:
  row[0] = random.random()

# output to disc
output_table = "output_table"
t0 = time.clock()
arcpy.da.NumPyArrayToTable(tab, output_table)
print time.clock() - t0

arcpy.management.Delete("in_memory")

I used this script on a table with 150K rows and 4 integer columns (~5MB filesize) stored on a SSD. Reading a binary file of this size takes less than 1 second in C++ or Python. However, in ArcPy CopyRows(...) takes around 5 minutes (measured with Python's time.clock()). For the output, the performance is similarly bad. I am quite puzzled about this.

I already tried other alternatives to the input (e.g. using UpdateCursor and da.UpdateCursor but I could not achieve similar times).
Furthermore, I read about running the script externally, but this makes things worse. The script above is as fast as I can get.

Am I doing something terribly wrong here? How come file-IO is so slow in ArcPy? Is there any solution to this or do I have to write a convenient IO by myself?

Edit: As hinted by @Jason Scheirer, the result of the internal computations can be written directly to disc, no need to go via internal memory again.

Update: I had the problem checked on a alternative environment by a client (organization with >50 Arc* users). In their native Windows environment the performance was equally bad. So I assumed the problem was not related to my system. However, it turned out that on my Windows inside a virtual machine the problem is solved by keeping the data on a local hard drive, dedicated to the vm. I never experienced performance issues with other applications accessing the parent (or a remote) file system. So I was misguided, sorry for that! I am curious what causes the bad performance on the client organization's environment... Will adapt the title and description ASAP.


texts: Extracting data from multiple raster layers in QGIS

texts:
I have a raster layer containing climatic monthly data:

rainfall
maximum and minimum temperatures

There is one layer for each variable and month.

I want to extract the average and maximum values for all of these layers into an Excel sheet.

Is there a way to do this without having to manually extract each value for every raster layer?


texts: How to get latitude, longitude and altitude database from a GeoTIFF file?

texts: I am developing an AR iOS application. When the iPhone gets the current location, it has to calculate locally the visible point. I want to read the tiff files in my application, and then get the coordinates and altitudes from it.(locally, offline) I am looking for an open source C++ code, or something like this. Please help


texts:
I have to convert from vector to raster so I used the Raster/conversion/rasterize tool. On my first try I used the Width (1447) and Height (2538) from the reference raster layer (1985clasification) on the raster size in pixel section in the rasterize window but the output layer differs on extent and pixel size. I tried a second time but changing the raster resolution in map.. based on the same reference layer and the output is the the same in extent from the first result but quite different in the way is rasterized.
I also tried to manually edit the gdal rasterize command using the -te parameter to see if I was able to match the extent but the output was completely different.

I followed the suggestion in a previous question but the extent of output layer keeps the same and it seems to lose information.

Via GRASS something is not working with my mapset so the vector is not listed and I decided to decline.

Am I using the wrong parameters ?  Do I need to include something else?

UPDATE

I managed to set the pixel size and extent. First the rasterize tool, then the raster calculator. Nevertheless the output (at least visually) is different. I assumed that the smaller raster (the one transformed from vector) is because the original does not have information beyond its own limits where the bigger raster does. How is possible for this two rasters to have same values and look like that?

2 UPDATE

Checking the images outside Qgis seem to have the extent it should have.

I was not able to see the full extent because those are null values? or because the way it renders with layers beneath them?

Now there is a white frame on both images (on the left). Why the process of resampling end up with this? is it possible to remove it?


texts: How to read geopackage tiles using Geotools

texts:
I want to read tiles from a geopackage file using Geotools (v. 17.0). I found the class GeoPackageReaderTest.java on GitHub and modified the example to read tiles from my file:

// try to read geopackage tiles (modified from https://github.com/geotools/geotools/blob/06e864481205da8af17a7e5083c7894d484785fa/modules/plugin/geopkg/src/test/java/org/geotools/geopkg/mosaic/GeoPackageReaderTest.java)
GeoPackage geoPackage = new GeoPackage(inputFile); // inputfile: http://ngageoint.github.io/GeoPackage/examples/rivers.gpkg
GeoPackageReader reader = new GeoPackageReader(inputFile, null);
GeneralParameterValue[] parameters = new GeneralParameterValue[1];
TileEntry tileEntry = geoPackage.tiles().get(0);
ReferencedEnvelope referencedEnvelope = tileEntry.getBounds();
Rectangle rectangle = new Rectangle((int) referencedEnvelope.getWidth(), (int) referencedEnvelope.getHeight());
GridEnvelope2D gridEnvelope = new GridEnvelope2D(rectangle);
GridGeometry2D gridGeometry = new GridGeometry2D(gridEnvelope, referencedEnvelope);
parameters[0] = new Parameter<GridGeometry2D>(AbstractGridFormat.READ_GRIDGEOMETRY2D, gridGeometry);
String tableName = tileEntry.getTableName(); // "rivers_tiles"
GridCoverage2D gridCoverage = reader.read(tableName, parameters); // throws java.lang.IllegalArgumentException: Envelope must be at least two-dimensional and non-empty.
RenderedImage img = gridCoverage.getRenderedImage();

In the line reader.read(...) I get an IllegalArgumentException:

  java.lang.IllegalArgumentException: Envelope must be at least
  two-dimensional and non-empty.

If I am printing the width and height of my envelope I get this result:

System.out.println("2 dimensions of gridEnvelope - width: " + gridEnvelope.getWidth() + " | height:" + gridEnvelope.getHeight());
// -> 2 dimensions of gridEnvelope - width: 4.0075016E7 | height:4.0075016E7

Any suggestions what is wrong with the code?


texts: Euclidean distance and connecting lines between polygons using ArcMap?

texts:
I have a large set of polygons. I want to

(1) calculate pairwise distances between polygons (so one polygon should be connected to several others). I can do this using the generate table in ArcMap.

and then

(2) I want a lines connecting these polygons (again one to many/many to one) that can link to the  distance file.

The only way I can think of doing this is through linkage mapper, but I dont want all the least-cost etc analysis, which is very time consuming.

Is there a simpler way to do this using ArcMap 10.5?


texts: Plotting points of restaurant locations onto NYC Taxi Data Shapefile Plot

texts:
I'm trying to plot the locations of restaurants in NYC onto the plot of NYC boroughs that I made from the shapefile provided by the NYC Taxi Data website using Rstudio. I don't understand why it is plotting the points of the restaurant locations wrong.

aoi_boundary = the shapefile provided by the NYC Taxi data website. (https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page)
rest_list = the list of restaurants with columns with the longitude and latitude coordinates where the restaurant is located

This is what I've attempted along with the attached result:

aoi_boundary_HARV <- st_read("/Users/rodrigotrujillo/Documents/College/Fall 2019/STAT 405/Project/taxi_zones/taxi_zones.shp")

ggplot(data = aoi_boundary_HARV) + theme_bw() +
  geom_sf(size = .2, color = "black", fill = NA) +
  ggtitle("Yearly Pickup Location Density Map") +
  coord_sf() +
  geom_point(data = rest_list,
             aes(x = Longitude, y = Latitude, color = GRADE,
                 size = 2))


texts: Calculating Polygon Width at Increments using ArcGIS for Desktop?

texts:
I am trying to measure the width of a polygon every km. I have a center line and a polygon, but am unsure what to do next.

In Estimating river Channel Width using Python/ArcGIS/MATLAB/R? I have read about methods which use create perpendicular transects and then using the intersect tool to measure the distance between the two sides of the polygon.

This seems like a reasonable solution, and I have tried it, but for some reason the "Perpendicular Transects" tool is creating lines which do not appear to be perpendicular to me.

I am using ArcGIS 10.2 for Desktop.


texts:
I'm planning to use the NetCDF file format for storing my time-series observation data on a server. I would like to use the "NetCDF-CF" conventions. Before I start writing any scripts for importing data from my SQL Server database to NetCDF, I'd like to make an estimate of the disk space requirements for the NetCDF file.

My data looks like this

1000 stations (each station has ID, Latitude, Longitude, Elevation)
10 quantities (for example, air temperature, air pressure,
precipitation..)
10 years of data with hourly time resolution, that is 10 years*365 days*24 hours =
87600 hours per station and quantity
the observation data value fits in 4-byte integer (if it's a decimal
number then I can scale the value to integer.) So I need 10 years*365 days*24 hours*4 bytes = 350400 bytes per station and quantity.
I presume that I also need to store the "time" information somewhere
in a suitable format. Usually 8 bytes is enough to store a time
value, so I need 10 years*365 days*24 hours*8 bytes = 700800 bytes for the time column.

Now calculating the NetCDF file size:

( [1000 stations] * [10 quantities] * [87600 hours] * [4 bytes per value]) + [700800 bytes for time] = 3504700800 bytes = 3.264 Gigabytes total

Is this roughly correct?

Or does NetCDF have any significant "storage overhead" like most relational databases do (for example in Postgres there is a 23 byte overhead per data row or in MSSQL it's 7 extra bytes per data row)?


texts: Deleted polygons seem to revert to older versions

texts:
In QGIS 2.12 carrying out the following steps leads to unexpected results. I wonder if anyone can shed some light on this.

Copy a layer from a previous project to a new project using Windows Explorer copy and paste of all layer files (dbf, prj, qpj, shp, shx).

Delete the polygons in the new project, make some new polygons in different locations and save.

Paste some data into the .dbf file of the new layer (being careful not to re-arrange the order of the polygons)

Returning to QGIS, the polygons have all reverted to their previous locations, this despite them being deleted, re-drawn and saved.
This is not consistently replicable, but I've not been able to work out what it is about the projects that do this.
Any ideas where QGIS is storing the locations of the old polygons so that I can clear this 'cache', or why it's doing this at all?
Edit
I've gone through a simplified set of steps to clarify the issue. Copy an entire project, including it's layers to another folder. Open one of the layers and delete the polygons. Save the changes, close the project. Open that layer's .dbf file, make one small change (I added the letter 'a' to a text column), save. Re-open the project - all the polygons have reverted to their previous locations before being deleted (even though I saved the deletion). How does QGIS know where the polygons used to be?
Edit 2
Sorry for the long edit, but this is really weird - Open the .dbf in OpenOffice Calc it has just the new polygons, open it in OpenOffice Base (which deals with dBase files) and it has just the new polygons. Open it in WPS (Excel equivalent), and all the old polygons are still there despite being deleted and saved twice!

Just to save anyone the bother, I've written this one off as a bug in QGIS, and so not really suitable for this forum.
The final proof - I loaded the offending layer in a friend's ArcGIS, all the polygons are still there, despite being deleted in QGIS twice.


texts: Assigning CRS to Lat/Lon in R giving error: "Geographical CRS given to non-conformant data"

texts:
I have a list of suburbs in Queensland, Australia along with their lat/longs.

However, when I try to read in the file in R and assign a CRS, R is giving the following error message:

  Error in proj4string<-(*tmp*, value = ) :
   Geographical CRS given to non-conformant data: 153.53332931

Here are relevant lines of the code,

# Loading suburb list
qs <- read.csv(file = '../data/Australian_Post_Codes_Lat_Lon.csv',
               header = TRUE,
               stringsAsFactors = FALSE)

suburb <- qs %>%
  filter(state == 'QLD') %>%
  select(postcode,
         suburb,
         lat,
         lon)

coordinates(suburb) <- c('lat', 'lon')
# Assigning proper Coordinate Reference System (CRS) to Lat/Long.
proj4string(suburb) <- CRS("+init=epsg:4326")

I am using the csv file provided here


texts: dragging while holding ctrl (control) key doesn't work in OpenLayers3?

texts:
When I hold down the  Ctrl (control) key on my keyboard and attempt to pan any OL3 map with my mouse, the pan doesn't happen.

Because of the particular needs of my application, I need panning to work when the ctrl key is held down, just as it does when Ctrl is not held down.

I couldn't see anywhere in the OL3 documentation about ctrl-drag being treated specially, so I don't know where to look to disable this behavior.

Any ideas?


texts: Stopping edit session at start of ArcPy script?

texts:
I have a simple script that I use to select features from a ArcSDE database and export them to a personal.gdb.

Currently, I have to manually save and stop my edit session in the MXD so that the script can get a schema lock.

I have been trying to use arcpy.da.Editor at the start of the script to automatically end the edit session, but I have not been able to achieve that thus far.

Is this the correct way to end an edit session in the MXD I am running the script in? I would also like to be able to restart my edit session after the script runs.

import arcpy
import os
from arcpy import env
env.overwriteOutput = True

edit = arcpy.da.Editor
edit.startEditing(False, True)
print 'Stopping Edit'

# Delete Feature class from file Geodatabase
arcpy.Delete_management('C:\GIS\BlairMailer\Mailer.mdb\AddressPoint')
print 'Feature Class Deleted'

# Select by attribute
arcpy.SelectLayerByAttribute_management ( "AddressPoint", "NEW_SELECTION", "   [DISTRICT] = 'Y'" or "[DISTRICT] = 'y'")

# Export selected feature to Geodatabase
arcpy.FeatureClassToGeodatabase_conversion("AddressPoint",'C:\GIS\BlairMailer\Mailer.mdb')
print 'feature Class Exported'

# Open Mailer
os.startfile(r'C:\GIS\BlairMailer\Mailer.docx')

This is the returned error when I try to run the script with an active edit session.

Runtime error
Traceback (most recent call last):
File "<string>", line 7, in <module>
TypeError: descriptor 'startEditing' requires a 'Workspace Editor' object   but received a 'bool'


texts: Export attribute table to Excel readable format

texts:
I want to save some attribute tables in a file format that Excel can read. Currently I'm using the tool Table to dBASE in Modelbuilder, but when I open the files in Excel it shows some errors in characters like "ç", "à", "õ". For example, what should appear like "Vila Chã" actually appear like "Vila ChÒ"

I need a format that accepts this type of characters correctly. How can I work around this situation? The solution should be a geoprocessing tool that can be used in a model, like "Table to Excel"

I've installed the tool "Table to Excel" and this is the error I'm getting:

There's some problem with encoding, but how can I solve this?


texts: Summing multiple routes using ArcGIS Network Analyst

texts:
I have been working with Network Analyst to try and solve the following problem, but the only solution I have found is tedious for more than about 4-5 origins and destinations.
I have a set of about 20 origins and 100 destinations connected via a network.  Goods should move from specified origins to specified destinations via a least cost (shortest distance) route along the network.  So, there are 2000 total routes. Each routes has a quantity of goods associated with it, for instance, Origin A sends 1000 kg of goods to Destination 1, 500 kg of goods to Destination 2, etc.
First, is there a way to process all 2000 routes at the same time?
Second, is there a way to sum the total amount of goods and total number of shipments that move along each part of the network?
For example, if 20 shipments carrying 15000 (total) kg of goods move through a high traffic area, can I weight that part of the resulting routes file to be thicker/heavier than a part of the network that receives 1 shipment carrying 100 kg of goods.  Therefore, congestion areas could be identified.  Ideally, I would like to be able to sum both the number of shipments and total goods transported along each part of the network.
I suspect I need to use ModelBuilder, but maybe there is a solution I am overlooking.


texts: Difference between GEOSERVER.war and GEOSERVER.exe?

texts:
Can someone comment on the performance/usage differences and advantages between deploying GeoServer from a web archive or installing GeoServer with the Windows installer?

I currently have GeoServer deployed in Apache Tomcat is this necessary to install another (GeoServer with the Windows installer)?


texts: Unexpected behaviour using ST_Equals() with QGIS

texts:
Let's say there is a table 'polygon_a' and a table 'polygon_b'. Each table contains three polygon features. The polygon features are spatially equal.

CREATE TABLE polygon_a (
gid serial NOT NULL,
geom geometry(polygon, your_SRID),
CONSTRAINT polygon_a_pkey PRIMARY KEY (gid)
);

CREATE TABLE polygon_b (
gid serial NOT NULL,
geom geometry(polygon, your_SRID),
CONSTRAINT polygon_b_pkey PRIMARY KEY (gid)
);

Now I've created a view 'equals_false' to select all features from table 'polygon_a' where ST_Equals(polygon_a.geom, polygon_b.geom) is false. After using the QGIS 'Split features' tool to split a feature from table 'polygon_a' the view 'equals_false' returns two results.

CREATE VIEW equals_false AS SELECT
polygon_a.gid,
polygon_a.geom
FROM
polygon_a LEFT JOIN polygon_b
ON ST_Equals(polygon_a.geom, polygon_b.geom)
WHERE polygon_b.gid IS NULL;

But after using the 'Merge selected features' tool to merge the splitted features 'equals false' still returns a result. What's the reason of this behaviour?


texts: How to select features that intersect but don't touch (using ArcObjects)?

texts:
So how do I do that? ArcMap's Select By Location selects also all features that touch, but I only need the ones that actually intersect the area.

I currently use esriSpatialRelEnum.esriSpatialRelIntersects and there doesn't seem to be an option for what I need. esriSpatialRelOverlaps isn't an option because it doesn't cover the case when the geometries are excatly the same size and over each other. esriSpatialRelCrosses also doesn't work because it doesn't support Area/Area intersection.

So is this possible at all?


texts:
i am using qgis 2.8.2 and qt4.8.5
i have a problem when i call zoomIn or ZoomOut, the raster on canvas is reduce or zoomed but :
- the raster size change
- it disaprar
- it appear again but 3 to 5 seconds later

i don't understand why there is this delay, it seems that there is 2 rendering of the mapcanvas

here is my code :

void CMyQgsS57::SlotzoomIn()
{
    ((CMyQgsMapCanvas*)mpMapCanvas)->zoomIn();
}

void CMyQgsS57::zoomIn()    //+
{
    if (((CMyQgsMapCanvas*)mpMapCanvas)->isDrawing() == true) {
      QTimer::singleShot(5, this, SLOT(SlotzoomIn()));
      }
}

void CMyQgsS57::SlotzoomOut()
{
    ((CMyQgsMapCanvas*)mpMapCanvas)->zoomOut();
}

void CMyQgsS57::zoomOut()   //-
{
    if (((CMyQgsMapCanvas*)mpMapCanvas)->isDrawing() == true){
        QTimer::singleShot(5, this, SLOT(SlotzoomOut()));
    }
 }

if you have any idea please let me know


texts:
I have been working on a mapproxy WMS/WMTS service for imagery. Everything has gone relatively smooth, except for when loading the WMTS service in ArcMap. Certain tiles are blurry and never adjust resolution when changing scale. However, the WMTS service works fine in QGIS and even ArcGIS Online.

What could I be missing in my configuration? Or is there another issue I am not aware of?

Link to the demo page: http://162.222.183.191:8080/demo/

EDIT:
WMTS GetCapabilities link: http://162.222.183.191:8080/wmts/1.0.0/WMTSCapabilities.xml

    services:
      demo:
      tms:
        use_grid_names: true
        origin: 'nw'
      wmts:
        restful: true
        restful_template: '/{Layer}/{TileMatrixSet}/{TileMatrix}/{TileCol}/{TileRow}.{Format}'
        kvp: true
        srs: ['EPSG:4326', 'EPSG:900913', 'EPSG:3857']
       wms:
         md:
           title: Mapproxy_Test
           abstract: Mapproxy_Test
         srs: ['EPSG:4326', 'EPSG:900913', 'EPSG:3857']
         image_formats: ['image/png']
         max_output_pixels: [4000, 4000]
     layers:
       - name: Mapproxy_Test_RGB
         title: Mapproxy_Test_RGB
         sources: [rgb_cache]
     caches:
       rgb_cache:
         grids: [webmercator]
         sources: [rgb]
         format: image/png
         request_format: image/png
         disable_storage: false
         concurrent_tile_creators: 5
     sources:
       rgb:
         type: tile
         grid: webmercator
         url: https://storage.googleapis.com/mapproxy_example/terrebonne-parish-imagery-pyramid-tiles/%(z)s/%(x)s/%(y)s
         transparent: true
         on_error:
           other:
             response: transparent
             cache: False
         coverage:
           srs: 'EPSG:3857'
           bbox: [-10221164.113,3320585.742,-10014524.117,3500502.768]
         http:
           ssl_no_cert_checks: true
     grids:
       webmercator:
         base: GLOBAL_WEBMERCATOR


texts: Add image in popup from json Url or image folder?

texts:
Added individual url and image for 21 places in json file.
Am trying to add image in my popup. Image shows up but broken in the popup. Not sure where in the code the error is??
(Same error for every popup, I rather load image from url, less space in project, but more control over pic if internal.)
My code:
var markerStyle = {
    radius: 3,
    fillColor: "#E4007C",
    color: "#E2062C",
    weight: 1,
    opacity: 0.9,
    fillOpacity: 0.7
};

var city = L.geoJson(cities, {
pointToLayer: function (features, latlng) {
    return L.circleMarker(latlng, markerStyle)
            .bindPopup(""
            + "
 "+features.properties.city + "
"
            + "County: " + features.properties.admin_name + "
"
            + "Population: " + features.properties.population + "
"
            + " GET HERE ");
    }
}).addTo(map);


texts: How am I supposed to use IWkb.ExportToWkb in C#?

texts:
ArcObjects provides an interface named IWkb for obtaining a WKB representation of a geometry. Generally, you can cast a geometry object to it:

// WKB representation of POLYGON((0 0,1 0,1 1,0 1,0 0))
byte[] wkb = Convert.FromBase64String("AQMAAAABAAAABQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADwPwAAAAAAAAAAAAAAAAAA8D8AAAAAAADwPwAAAAAAAAAAAAAAAAAA8D8AAAAAAAAAAAAAAAAAAAAA");

IGeometry g;
int bytesRead;
IGeometryFactory factory = new GeometryEnvironmentClass(wkb, out g, out bytesRead);

System.Diagnostics.Debug.WriteLine(g);
// Output: System.__ComObject
System.Diagnostics.Debug.WriteLine((IWkb)g);
// Output: System.__ComObject
System.Diagnostics.Debug.WriteLine("Yes, the cast is valid");
// Just verifies we got past the cast without error.

However, the signature of IWkb.ExportToWkb doesn't seem to make sense. The docs show takes two parameters:

ref int byteCountInOut
ref Byte& byteBuffer

There are three things that don't make any sense:

I know ref means that it takes an incoming value and might modify that value to also use for output, but sending input for these variables doesn't make sense; I'm trying to get output from them.
byteBuffer is of type Byte& according to the docs, but that trailing & doesn't appear to be valid C# syntax.
When I look at the generated metadata for the method in Visual Studio, it shows that byteBuffer's declaration is actually out byte byteBuffer. It being an out parameter makes more sense than it being ref, but the fact it's byte makes even less sense. A WKB consists of many bytes, so it needs to be some data type that allows for that.

How is this method supposed to be used? Do I need to use some kind of loop to get one byte at a time; if so, how? Is there a better interface for getting the WKB form of a geometry?


texts: Open source location allocation modelling

texts:
I'm looking for open source tools for location allocation modelling (eg planning public health facilities).  Preferably in a GIS environment.  I know GRASS somewhat, but v.net.alloc doesn't seem to take account of demand.

I don't know if SAGA or Sextante have something?

Thanks!


texts: How to hand over polygon attribute to new point while digitizing?

texts:
My task is to work with bird observation data in an area, that is divided into zones as different polygons (A1, A2, etc.).

I have the hand-written lists of observations and want to bring all entries into a database. My idea was to create a point feature somewhere in the corresponding zone for each line, and then enter all data that is provided (Species, number, begin of observation, end of observation, etc ...) with a custom form.

I also need to have the polygon ID in the data, which I of course could enter manually, but I would rather prefer to learn how to include it automagically by querying to position of the point over the polygon and extract the polygon ID somehow.

As I want to work with a spatialite database, i thought about creating a trigger that queries the polygon the new point is created above and include the ID in my point data.

My problem: I have no idea how to start and can't find examples on that.

Anyone some ideas?

Thanx Bernd


texts: Modelling per-segment attributes in GML LineString

texts:
This is not a software-specific question, but rather a question about how to represent some concepts in a GML LineString.

It has to be GML, for policy reasons.

Consider a GML LineString consisting of contiguous segments. I'm trying to represent attributes that are valid for a segment, but not necessarily valid for  the whole LineString.

As an example, it might be a car route along several roads. Each road on the route could have some different properties (surface, minimum and maximum speeds, width, surface roughness). My modelling scenario is quite different to the car route, so I'm not looking for something that is specific to that domain (e.g. the OSM model isn't going to work for this - it is not GML for a start).

Conceptually, the LineString could have some default attributes, which are overridden by per-segment attributes. In the car route example, the route might have a default speed, but some segments are higher or lower speed.

The concepts I've explored so far:

Using a LineString, and adding numbered attributes / elements.
Using a MultiLineString or collection, where each LineString is a single segment, attaching the attributes / elements to a LineString.
Using Points to implicitly represent the LineString, such that the attributes / elements on the Point apply to the next Point.

None of those seems very satisfactory (loss of explicit representation, possible gaps in the collection cases, non-obvious intent).

Is there published guidance on performing this kind of modelling (e.g. from OGC or ISO), or failing that, anything that is likely to be better supported in common tools / GML profiles than what I've already identified?


texts: Accessing ArcSDE delta tables through Python?

texts:
I would like to access delta tables through ArcPy for a QC script that I am improving. Read-only. How could I output a Python dictionary that contains all the rows that were added/deleted/modified between a named version and DEFAULT for a particular feature class?

We are currently on an Oracle client, but we will be switching to SQL Server, if there are substantial differences between the two.


texts: Layer from PostGIS with array of integers fields?

texts:
I have GeoServer 2.6 and I published a vector layer (point) which some of the field attributes are integer's array.

GeoServer doesn't recognize those fields in wfs call or even in the Feature Type details table.

Is there any possibility to give them visibility?

I need to symbolize that points taking in account those fields.


texts: I am trying to create a choropleth map of ice concentration over Lake Erie. I have data for all the lakes (I can clip it all at the end), which are in a point shapefile. I have laid a grid over the area and I am able to create a union between the two, but the attribute table pulls up NULL for each of the ice data attributes added to the new grid table. Data can be found at: http://www.glerl.noaa.gov/data/pgs/glice/glice.html (I used 20080211) I am just not sure how to assign a value to each grid cell with the given data type.


texts: Area of parcel on Autocad is deferent from ArcGIS and we need to rectify ArcGIS Area Calculations

texts:
I have a problem with the built in Area calculation field in ArcGIS and in Autocad map 3d area for the same feature. The issue happened specifically when using Arcs (Curves) .

The scenario is we have drawn using Autocad map 3d then we redraw it on ArcGIS in a different way, sometimes using trace and sometimes import or manual drawing. The issue is these different methods result in different area calculations between the ArcMap calculation field and Autocad query.

Any Ideas?


texts: Finding available options when adding layer to the WebMap() model through ArcGIS Python API?

texts:
I'm trying to find available options I can use when adding a layer to WebMap() in the ArcGIS API. The esri examples don't expand on all the available options and when looking at the API reference located here https://developers.arcgis.com/python/api-reference/arcgis.mapping.toc.html#arcgis.mapping.WebMap.add_layer it doesn't expand on all the available options either.
Is there a location where I can find all the available options?


texts: Match attributes (with spatial join) weighted by length of line feature in polygon in ArcGIS

texts:
I have line features (streets) and polygons (air quality raster). I need to transfer the air quality attributes into the lines. However, some lines are part of two or more polygons (e.g. the selected feature below).

I am looking for a way to transfer the attributes of the polygons onto the lines weighted by the length of the part of a line feature within each polygon.

The tool Spatial join only allows for a simple average in its merge rules (mean) giving equal weights for both parts of the selected line in the picture above.

An alternative idea was to split line features at intersections with polygon boundaries and have the air quality attributes assigned to the two or more new features. In a second step, I could then unsplit those split features using an ID. Lastly, I could use the field calculator to calculate a new weighted average air quality field, like:
Air Quality = (Shape_Length A * Air Quality A) + (Shape_Length B * Air Quality B)/(Shape_Length A + Shape_Length B)

I am using ArcGIS Desktop 10.5.


texts: Adding CSV with state plane coordinates instead of decimal degrees to ArcGIS online?

texts:
My goal is to use this method of using google sheets with ArcGIS online to update my web map automatically with new points which are added to the google sheet csv.

My problem is according to this document, coordinates have to be in decimal degrees. The points I want to bring in are surveyed in state plane northing and easting. It would defeat the purpose of updating it automatically if I have to create a shapefile out of the csv first in order to get decimal degrees fields.


texts: Using conditional statements in ArcGIS Field Calculator?

texts:
I have two fields, "P_Area" and "PC_Area", whose sizes I'd like to compare in a third field, "Areas_Equal". The "Areas_Equal" field should contain a "Yes" if the other two are equal and a "No" if they aren't.

How can I do this with Field Calculator and its Python Parser?


texts: How to extract height of building and roof structure from nDSM data

texts:
I am working with QGIS; SAGA; OrfeoToolbox and I try to create LOD1 and LOD2 building models.
I need help to recreate LOD1 (or LOD2) building models from DSM / LIDAR data.
I have a footprints of buildings and both DTM, DSM and raw LIDAR Data.
My LOD1 workflow and problems
If building is homogenous then I know that I can use nDSM and ZonalStatistics with Max value to achieve LOD1. But what if in one building footprint there are many building parts with different heights? I attach SHP footprints and nDSM data for my buildings. You can see clear differences in building parts in both left and right building. I have to proceed 100 000 buildings that why I search for quite autonomous solution.

My Idea - Is there any way to create height break-lines in nDSM data - and then to run polygonize function and ZonalStatistics?
Any other suggestion?
And now about LOD2
I decided to create dense diamond fishnet, intersect with building footprints, and drape it on nDSM with 3D-multipolygon. But results don't look good enough. And I can't find way to simplify roof structures with similar angle. I wonder if there is any access to typ-of RANSAC generalisation function to simplify that areas?
What I have now:
Any other suggestions and tips? how to achieve as accurate as possible LOD2 models inside QGIS?
Will be glad for any suggestions and tips.


texts: How to calculate values in a 3*3 matrix using raster calculator?

texts: How to manipulate a raster values in a 3 X 3 matrix.. using raster calculator or any other method e.g. to find out how many cells are there in a 3 X 3 matrix having cell values less than the center cell.. and sum them up...


texts: How can I prevent the Legend from covering each Atlas Feature?

texts: I'm generating an atlas with a set of polygons. I want the legend to be above the map, but to ensure that the legend doesn't actually cover the featured polygon. Is there a way to do this?


texts: ArcGIS Toolbox functions missing all parameters?

texts:
I recently installed ArcGIS 10.3.1 and now all the toolbox functions (that I have tried) open, but they have no parameter fields. For example, the "Add Field" tool comes up with the standard toolox dialog with the title bar, buttons at the bottom, and a functional help panel, but the rest of the dialog is entirely empty.

I've run a repair on the installation.

Has anyone else seen this? Any fixes?

Same thing after uninstalling all ArcGIS applications and reinstalling.
Other details: Windows 7, 64 bit, SP1, Intel CPU


texts:
I am looking for a Python code/function that can extend an arc linestring up to a pre-defined length while maintaining the same arc of the original line. The idea is to repeat the configuration of the last 2 segments of a linestring n times. My first inclination was to get the angle between these 2 last segments and apply affinity.rotate using always this angle while updating the origin in a while loop, something like:
from shapely.geometry import LineString, affinity
import numpy as np

n=100
last2SegmentsOfaLinestring=[(3,0),(3,5),(7,8)]
angle = alpha #angle between[(3,0),(3,5)] and [(3,5),(7,8)]
line = LineString([(3,5),(7,8)])
newX, newy = 7, 8
i=0
listOfTuples = []
while i < n:
    rotated = affinity.rotate(line, angle, origin=(newX, newy))
    # convert rotated to list of tuples to get updated (newX, newy) - this is the line centroid
    listOfTuples.append((newX, newy))
    i+=1

If I am thinking this correctly, I will need to convert the list of tuples to linestring, apply affinity.rotate, and convert back to a list of tuples to get the new origin. However, this seems a bit computational expensive and I wonder if there is an optimized function or a better way to do this.
Another way I tried is to apply a rotation to the vector and keep track of the origins of the next segment, adding the new, coordinates. This seems a bit more straightforward.
orignalVector = np.array([0,5])
rot = np.array([[np.cos(rotAngle), -np.sin(rotAngle)], [np.sin(rotAngle), np.cos(rotAngle)] ] )

xx = np.dot(rot, orignalVector)

but the question remains. Is there an optimized function that does all that?


texts: Find geometry of shortest paths of series of points in a shapefile

texts:
I have two layers. The first where there is a set of routes, and another layer where there is a series of points. I would like to find the shortest paths following the route layers.
I used OSM route on qgis but it is limited on only two points.
I am looking to automate OSM route by calling on my two layers.


texts: Fill color in polygons using attribute variable in GMT (Generic Mapping Tool)

texts:
I want to FILL colors of the polygon (converted from .shp ESRI shapefile) based on an attribute (here OBJECTID, Area) from lookup color (.cpt), but I am getting colored outline of each polygon. I attached my result (polygon.png) using below command.

psxy river_b.gmt $range $proj -X2i -Y3i -Ba1WESN -aZ=OBJECTID -Criver_b.cpt -K -P >! $outfile

Please suggest me which gmt command will work here.

I am attaching a example what kind of map I want to make in GMT.


texts: CartoDB displaying layers and infowindow?

texts:
I just started testing and using CartoDB, and I am dealing with this problem on a map that has been published already:

http://deix.is/ueberwachungsstatistik/

I have made minor corrections to the custom HTML infowindow code, and now

a) the infowindows in FF span vertically to the extreme.

b) all layers appear as "visible", even though only 2 of 4 layers should be - in the CartoDB interface, only 2 are visible, the rest is hidden. It also doesn't matter which one.

Maybe the two issues are connected? I have tried lots of refreshing, checking and repasting code, and really don't know a way out.


texts: I have a problem to display the legend in QGIS Web Client

texts:
I have installed QGis Map Server and QGis Web Client on Centos 6.4 64-bit operating system, everything works fine except the legend box which is filled with garbage data instead of text (it seems that some fonts are missing), please have a look of the screenshot attached.

OS:
CentOS release 6.4 (Final)

Kernel:
2.6.32-358.2.1.el6.x86_64 #1 SMP

Software packages:
qgis-mapserver-1.8.0-9_0.el6.x86_64;
qgis-1.8.0-9_0.el6.x86_64;
mod_fcgid-2.3.7-1.el6.x86_64;
fcgi-2.4.0-10.el6.x86_64;
httpd-2.2.15-28.el6.centos.x86_64;
postgresql92-libs-9.2.4-1PGDG.rhel6.x86_64;
postgresql92-9.2.4-1PGDG.rhel6.x86_64;
postgresql92-server-9.2.4-1PGDG.rhel6.x86_64;


texts: How to remove spaces in only the top row in excel?

texts: I want to clean up some data that's having trouble importing into Arc.  I want to know a way to clear spaces in only the top row in excel.  Sometimes there is one hidden before or after the row name and it causes Arc to import them incorrectly.  It will return an invalid SQL statement when it tries to draw in Arc.  I'm using ArcGIS 10.2.2.


texts: ESRI ArcGIS Runtime for .NET view non ESRI Vector Tiles

texts:
I am trying to view some vector tiles hosted outside of ArcGIS Online/ArcGIS Enterprise on my own tile server (Flat non indexed tiles).
I have been unable to add the Tile Server to the .NET Runtime application.
I think it's because of ESRI's implementation of the Stylesheet baked into the URL.

ESRI wants {level}/{Col}/{row}.pbf  instead of standard placeholders /{z}/{x}/{y}.pbf  (at least that's how you add tile server with the different placeholders)
https://tiles.arcgis.com/tiles/a6SBLWoOui928qXS/arcgis/rest/services/UAS/VectorTileServer
the stylesheet is found by default with ESRI https://tiles.arcgis.com/tiles/a6SBLWoOui928qXS/arcgis/rest/services/UAS/VectorTileServer/resources/styles

the URL for ESRI's basemap
https://basemaps.arcgis.com/arcgis/rest/services/World_Basemap_v2/VectorTileServer


texts: How to write natural logarithm in SNAP software

texts:
I want to compute the temperature of a band using the equation:

where L(6) is Band 6 of a Raster file and K1 and K2 are constant.
I am using Band Math to perform conversions. However I could not find a natural logarithm under the Functions in Band Math Expression Editor.

There are functions for log and log10. How do I go about this?


texts: Labeling without intersect geometries in QGIS?

texts:
I have a polygon geometry and I have some points that are its vertices.
QGIS 2.18 changes randomly the position of labels. When I choose 'outside of symbol' option, labels are showed outside of point symbol, but not outside of polygon.

I am doing many maps and I need labels to be showed like this.

How could I show labels of these points changing its position to be outside of the polygon automatically?

Is there any rule based solution or anything else?


texts: Save the result of raster calculator as an memory layer in pyqgis

texts:
I am using QGIS raster calc as below. Now I am saving it as a tiff. But is there a way that I can save it as a memory layer?

I want to refer to that resulting memory layer and polygonize it, and upload it to postgis server. Because I am building a plugin so I don't want user to specify a folder to save the tif.

layer = QgsRasterLayer(path_of_image, 'input')
name = layer.name()
resultinglayer = QgsRasterCalculatorEntry()
resultinglayer.ref = name + "@1"
resultinglayer.raster = layer
resultinglayer.bandNumber = 1
entries = [resultinglayer]
calc = QgsRasterCalculator('(%s@1)*0' % name,
                          'C:/tom/ras.tif',
                          'GTiff',
                          layer.extent(),
                          layer.width(),
                          layer.height(),
                          entries)
calc.processCalculation()


texts: Receive "ValueError: Number of label classes don't match number of class breaks" when trying to update raster symbology with arcpy

texts:
I'm working on an arcpy script that will automatically set up excluded values, class breaks values, and class break labels in a series of raster layers. A simplified version of part of my code is visible below:

mxd = arcpy.mapping.MapDocument(r'C:\Users\djh\Desktop\GIS Temp\test.mxd')
df = arcpy.mapping.ListDataFrames(mxd)[0]

symbology_lyr = (r'C:\Users\djh\Desktop\GIS Temp\test.lyr')
raster_name = (r'C:\Users\djh\Desktop\GIS Temp\test.gdb\test')

levels = [31.0, 32.0, 33.0, 34.0, 35.0]
class_break_labels = ['31.0 - 32.0', '32.0 - 33.0', '33.0 - 34.0', '> 34.0']

layers = arcpy.mapping.ListLayers(mxd, data_frame=df)

for lyr in layers:
    if lyr.isGroupLayer:
        for sublyr in lyr:
            if sublyr.name == raster_name:
                arcpy.ApplySymbologyFromLayer_management(sublyr, symbology_layer)
                if sublyr.symbologyType == "RASTER_CLASSIFIED":
                    sublyr.symbology.excludedValues = '0 - ' + str(levels[0])
                    sublyr.symbology.classBreakValues = levels[1:]
                    sublyr.symbology.classBreakLabels = class_break_labels

When I run up to the Apply Symbology From Layer line, it works fine and updates the symbology accordingly, but when I include the last 'if' loop, I get the following error on the last line of the above snippet:

  ValueError: Number of label classes don't match number of class breaks

I included a print statement to get the difference between the lengths of levels[1:] and class_break_labels (both lists), and for every iteration it was 0, so I don't see how the numbers could possibly not match.

I appropriated the last 'if' loop from this Esri help page.

Has anyone encountered this before?


texts: Creating new Geodatabase, FeatureClass, Table and adding Fields using ArcGIS Pro .NET SDK

texts:
I need to create a new Geodatabase and creating FeatureClass, Table in the Geodatabase. After that adding fields in FeatureClass and Table are also needed.
I have created a new Geodatabase using this code:
     private async Task<bool> ExecuteAddFileGDB(string fileGdbPath, string fileGdbName)
  {
        try
        {
            return await ArcGIS.Desktop.Framework.Threading.Tasks.QueuedTask.Run(() =>
            {
                var fGdbPath = fileGdbPath;
                var fGdbName = fileGdbName;
                var fGdbVersion = "Current";
                var parameters = Geoprocessing.MakeValueArray
                    (fGdbPath, fGdbName, fGdbVersion);
                var cts = new CancellationTokenSource();
                var results = Geoprocessing.ExecuteToolAsync("management.CreateFileGDB", parameters, null, cts.Token,
                    (eventName, o) =>
                    {
                        System.Diagnostics.Debug.WriteLine($@"GP event: {eventName}");
                    });
                return true;
            });
        }
        catch (Exception ex)
        {
            MessageBox.Show(ex.ToString());
            return false;
        }
   }

I have also code developed for creating FeatureClass, Table, and for adding Fields using Geoprocessing Tools.
As my ArcObjects experience, I think a geoprocessing operation is comparatively costly than an API call. So I am seeking a way to create a new geodatabase using ArcGIS Pro API instead of the Geoprocessing tool. Like ArcObject API http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//0001000004t8000000.
As well as for creating FeatureClass, Table, and for adding Fields.


texts:
I am trying to delete clouds (ideally also their shadows) from Rapid Eye images (5 bands; blue, green, red, red edge, near-infrared).

I have tried to set a threshold for the red band and mask accordingly, inspired by this question but it doesn't pick up the clouds, but only other areas (such as parts of a lake and vegetation and also a few of the clouds shadows). I have also tried to modify the algorithm: ee.Algorithms.Landsat.simpleCloudScore() without success (don't know if it would even be possible to modify landsat algorithms and apply to other images).

Is there any other convenient way to mask the clouds?


texts: Is it possible to detect human settlements by night light data?

texts:
I am working on a project in which human settlements are identified using geo-tagged Twitter data. Since this project is focused on developing countries (mostly African), I am afraid that many human settlements may not be identified because Twitter may not be used there much.

I have been thinking about using night light satellite imagery instead. My assumption is that if there is Twitter data, there is electricity and therefore there is likely to be streetlights or other light pollution. I would like to try this method and compare the results.

Question: Are there high-resolution night light datasets available that have the potential to produce accurate information about human settlements? Are there perhaps some drawbacks or problems with this approach that would preclude this method from success?


texts:
I am creating a PNG that is referenced to the OSGB36 datum (1 pixel per 50m easting/northing that is colored according to the result of some calculation code)

I want to display this on Google Earth using KML. But obviously I need to transform the image before I do this.

I can manually use gdal_translate to create a TIF file then gdalwarp to convert the image but I want to achieve this in code (C++) not via external applications.

I looked at the GDAL API but having no joy doing what I can easily do with the aforementioned GDAL utilities.

http://wtp2.appspot.com/warpdemo.htm seems to doing something similar to what I want but having looked at the JavaScript I don't really understand how the affine transformation is being constructed.

Any pointers on how to achieve this would be much appreciated.


texts: Geoserver arabic labels and attributes doesn't show in right encoding character set

texts:
I set my shapefiles UTF-8 encoding which supposed to display Arabic fine, but it doesn't work in QGIS unless I change the encoding to Windows-1256. QGIS reads the arabic well with this encoding.

But when I publish to Geoserver using the plugin or by uploading the shapefile and specify the exact same encoding or connected via Postgre it doesn't work in layer preview or GeoExplorer or when I add the WMS into ArcGIS, for example.

Does anyone have any recommendations for eastern languages encoding especially arabic? Why is this happening while it should work fine?


texts: PostGIS Column has Z dimension but geometry does not

texts:
I am able to create a geomtry column for 3d multiline using this:

    CREATE TABLE trajectories
(
  gid serial NOT NULL,
  vehid character varying(10),
  state integer,
  geom geometry(MultiLineStringZ),
  CONSTRAINT trajectories_pkey PRIMARY KEY (gid)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE trajectories
  OWNER TO postgres;

But when I tried to alter an existing table by modifying an existing geom column or adding a new one, I am getting he following error:

    ERROR:  Column has Z dimension but geometry does not

  ****** Error ******

  ERROR: Column has Z dimension but geometry does not SQL state: 22023
  ERROR:  Column has Z dimension but geometry does not

ALTER TABLE cl ALTER COLUMN geom TYPE geometry(MultiLineStringZ)

How do I modify and existing geometry column to geometry(MultiLineStringZ)?
How do I add new geometry(MultiLineStringZ) column to an existing table?


texts: How to apply a Cloud Mask in a collection?

texts:
I am trying to apply a mask cloud in a SAVI collection (LANDSAT) to get a chart with time variations between pixels.

var region = table3

    Map.addLayer(table3)

    var col = ee.ImageCollection('LANDSAT/LE7_SR')
        .filterDate('2004-01-01', '2011-12-31')
        .filterBounds(table3)
        .filter(ee.Filter.eq('WRS_PATH', 1))
        .filter(ee.Filter.eq('WRS_ROW', 72));

        print(col)

    // Mask clouds
    var col_noclouds = col.map(function(img) {
                       var mask =
    img.select(['cfmask']).neq(4).neq(1).neq(2).neq(3)
                       return img.updateMask(mask)
                       })

    // Median image
    var medians = ee.Image(col_noclouds.median())

    var clipper = medians.clip(table3)

    // A function to compute Soil Adjusted Vegetation Index.
    var SAVI = function(image) {
      return clipper.expression(
          '(1 + L) * float(nir - red)/ (nir + red + L)',
          {
            'nir': image.select('B5'),
            'red': image.select('B4'),
            'L': 0.5
          }).set('system:time_start', image.get('system:time_start'));
    };

    // Shared visualization parameters.
    var vis = {
      min: 0,
      max: 1,
      palette: [
          'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
          '74A901', '66A000', '529400', '3E8601', '207401', '056201',
          '004C00', '023B01', '012E01', '011D01', '011301'
      ]
    };

    var image = ee.ImageCollection(col)
       .map(SAVI);

    print(image)

    Map.addLayer(image,vis);

    var TS5 = Chart.image.seriesByRegion(image, region,
    ee.Reducer.mean(),'constant', 500, 'system:time_start').setOptions({
              title: 'NDVI Long-Term Time Series',
              vAxis: {title: 'constant'},
    });
    print(TS5);

I can make the chart for "image" (collection with SAVI function) but I do not know how to insert the cloud mask.

Any recommendation?


texts: In QGIS is there any way I can connect to and list tables for only one specific schema (among many schemas and not the same schema as my username) in oracle?  In the "Create New Oracle connection" box, can I use these two fields "Options" and "Workspace" to achieve this?


texts: calculating the population of destinations

texts: I would like to calculate the population of my destinations I created with OD Cost Matrix. I have got population data but I don't know how to find the population of my destinations scattered on my population data. Can you help me, please?


texts: Find square with latitude and longitude for given lat and lng and distance

texts: I have the coordinates of some point, lets say lat and lng,  I want to find the latitude and longitude of top left and bottom right corners of the square, whose sides are located k kilometers from the initial point.


texts:
I use QGis 2.4 over linux ubuntu 12.04
(a very standard situation :D)

Why does the list of installed plugins not appear in the program?


texts:
Can anyone tell what is wrong with this code?

import h5py, gdal, osr,gdalconst
from gdalconst import *

path3='D:/state_subsets'
directory=os.fsencode(path3)
for f in os.listdir(directory):
    filename = os.fsdecode(f)
print(filename)
g=gdal.Open(filename,GA_ReadOnly)
print(g)

the print(g) gives a none output.