texts: Max retries exceeded with URL in Nominatim with GeoPy

texts:
I'm trying to get a country name by lat lon (reverse encoding) using Python. I got to Nominatim with GeoPy which worked great two times but now I keep getting Max retries exceeded error.
My code:
from geopy.geocoders import Nominatim

def get_countries(lat,lon):
    coordinates = f'{lat} {lon}'
    locator = Nominatim(user_agent='myencoder', timeout=10)
    location = locator.reverse(coordinates,language='en')
    country = location.raw['address']['country_code'].upper()
    return [country,lat,lon]

get_countries(32.782023,35.478867)

In practice I would like to use this function for about 100-200 locations, but I get the error also for this one location. It worked the first time I run the code, but then I'm keep getting:
GeocoderUnavailable: HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Max retries exceeded with url: /reverse?lat=32.782023&lon=35.478867&format=json&accept-language=en&addressdetails=1 (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x000001F7AA229978>: Failed to establish a new connection: [WinError 10061] No connection could be made because the target machine actively refused it',))

I've tried changing the user_agent by writing my email or another name as suggested here, I also read the answers here and here, also other solutions. All resulted in this error.
What else can I do? Are there other reliable options using Python to get a country name (preferably with to letters e.g., ES, IT, US, etc.) by lat lon?


texts: Chosing a GeoCMS server collaboration tool

texts: I have seen a few things about different GeoCms. Researching these, none of them by default to actually seem to be a GeoCMS. Somehow when you hook these up to mapserver they become a GeoCMS. I do not quite understand how this exactly works and I am wondering if there are any good GeoCMS's that are good for digital collaboration (for example if you cannot speak with the user in real life) and something easy for the non-techy. I have looked at Wordpress, TikiWiki and Drupal but do not understand how they are GEOCm's


texts: Adding Raster to Weighted Overlay causes large parts of Output Raster to have NoData pixel values?

texts:
I've been building a large suitability analysis model in ModelBuilder over the past couple of months, it has a total of 11 rasters created and put into a weighted overlay. The model runs and an acceptable output is created when I have 10/11 Rasters used in the weighted overlay, I have one raster that for some reason causes the output raster to have Nodata Pixel values for large areas.
I can't give a more detailed description without typing a pretty lengthy write up so I hope this is enough, if not I will update the post when I have a chance.
Does anyone know why this might be happening?


texts:
I am  using the gimms package to automatically download NDVI data from 1981-2013. The code is provided here. The data is available here. However, after running it shows that data has been downloaded but nothing in the download folder.
What might be the problem?

## load 'gimms' package
library(gimms)

gimms_files <- updateInventory(sort = TRUE)

gimms_files[1:20]

## download entire gimms collection
downloadGimms(dsn = paste0(getwd(), "/data"))


texts: ArcGIS SDK for .Net (used in WPF): update map's zoom parameters after changing base layer

texts:
I'm having the following issue:
Suppose I have a Map object and 2 base layers (a world map and a street map). A world map doesn't let us view the map too closely. The street map allows us to view as close as we want (almost). So, they have different scale parameters. Now, if I load the world layer first, like this (pseudocode):

map.Layers.Insert(0, WorldLayer);

and then switch layers:

map.Layers.Remove("WorldLayer");
map.Layers.Insert(0, Streetlayer);

The map doesn't allow me to get closer than a WorldLayer allowed. SO, actually, I can't view a street view normally.

BUT, if I load a detailed map first:

map.Layers.Insert(0, Streetlayer);

and then switch layers:

map.Layers.Remove("Streetlayer");
map.Layers.Insert(0, WorldLayer);

In this case I can freely view trees near my home :) And even if I switch back, the WorldLayer would allow me to get closer than it allowed before.

I suppose there should be a way to update a map, to renew scale/zoom parameters for a new set of layers. Is anyone aware of a way to achieve that?


texts: Including use of Mask in Raster > Vector conversion in QGIS

texts:
The Polygonize (raster to vector) tool in QGIS Raster > Conversion menu allows for the use of a mask, however this option is not available in the Processing Modeler. What would be the best way to incorporate this functionality into the model?

For some context, I am converting a single-band raster to polygons using itself as a mask (values range from 0-1).


texts: It's possible to use gdal_calc.py with an entire path of rasters?

texts:
I try to use gdal_calc.py to scale and transform a time-series of MODIS images.
I've done that when i worked with a single image in the -A variable and a mask in -B, but i need to work with a directory of 350 .tif images and convert to ENVI .bin raster.
My problem is how to assign the -A variable to realize the operation with all the images and obtain 1 output file per image

To work with a single image i use:

gdal_calc.py -A image.tif -B MASK10000 --calc="A/B" --NoDataValue=0 --format=ENVI --outfile=image.bin

I use the linux bash, but maybe a python scrip in linux runs better


texts: Subdividing a GeoTIFF into 1x1 tiles of longitude and latitude?

texts: I have a big GeoTIFF file of the whole Slovakia territory in 25m resolution. It has around 17700 x 8000 pixels. I need to subdivide it into smaller files and every file must have 1° longitude and 1°latitude. I tried gdal srcwin but I don't know how many pixels is a 1° longitude and latitude and projwin didn't work for me. How do I find out how many pixels are 1° longitude and latitude?


texts: ArcGIS project raster produces incorrect vertical coordinates

texts:
I have a DEM in UTM (meters) / NAVD88 (meters) and need to transform it to State Plane (feet) / NAVD88 (feet). Unfortunately, I do not have the raster calculator in ArcGIS, so I have to use "Project Raster" and apply a vertical transformation. My xy coordinates are correct and the scale of my vertical coordinates is correct. Since they are both on the same datum (NAVD88) and only differ in scale, you'd think that the vertical transformation would be easy (just multiply the vertical values by ~3.28) and would work correctly. However, vertical coordinates shift upwards by ~ 5' from where they should every time I apply this transformation. Any idea why this happens?

Side note: I have resolved this problem in QGIS, but would like to try to figure out why I am having trouble in ArcGIS as I would prefer to perform all of my GIS operations in ArcGIS. The attached image shows the default transformation which is applied (unnecessarily/erroneously) during the operation.


texts: TypeError when entering info in leaflet-geosearch search bar

texts:
I am using leaflet-geosearch. When I begin typing anything into the search bar, nothing pops up, and the developer console reads:
SearchControl.ts:28 Uncaught (in promise) TypeError: Cannot read property 'label' of undefined
    at resultFormat (SearchControl.ts:28)
    at resultList.ts:38
    at Array.forEach (<anonymous>)
    at t.e.render (resultList.ts:35)
    at SearchControl.ts:345

Here is the code snippet:
  <link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css"/>
  <script src="https://unpkg.com/leaflet@1.7.1/dist/leaflet.js"></script>

  <!--initialize geosearch -->
  <link  href="https://unpkg.com/leaflet-geosearch@latest/assets/css/leaflet.css" rel="stylesheet" />
        <script src="https://unpkg.com/leaflet-geosearch@latest/dist/bundle.min.js"></script>

var GeoSearchControl = window.GeoSearch.GeoSearchControl;
    var OpenStreetMapProvider = window.GeoSearch.OpenStreetMapProvider;

var provider = new OpenStreetMapProvider();

var searchControl = new GeoSearchControl({
  provider: provider,
});

map.addControl(searchControl);


texts: Create polyline using MapBasic (unknown number of nodes)

texts:
I'm trying to create pline using MapBasic. My idea was to create simple tool with dialog box containing couple of EditBox, PopupMenu's to retrieve pline attributes from user. So far I have only making tools with same logic for point object and my code would look something like:

'************************************************************

Sub Main()

    Create ButtonPad AppName
        As toolbutton
        Icon 16
        DrawMode DM_CUSTOM_POINT
        Cursor MI_CURSOR_LRG_CROSSHAIR
        Calling SomePointSUB
    Show
    Fixed

End Sub

Sub SomePointSUB

    click_x = CommandInfo(CMD_INFO_X)
    click_y = CommandInfo(CMD_INFO_Y)

    Insert Into "POINT_TAB" ("attribute1", "attribute2"...., obj)
    values ((value1), (value2)....., CreatePoint(click_x, click_y))

End Sub

'************************************************************

So, this was straight forward, but creating polyline is not so (at least for me). All of the examples I found online had fixed number of nodes, like this (from MapBasic user guide):

'************************************************************

    Include "mapbasic.def"
        Type Point
        x As Float
        y As Float
        End Type

    Dim objcoord(5) As Point
    Dim numnodes, i As Integer, myobj As Object
        numnodes = 3
        set CoordSys Earth
    objcoord(1).x = -89.213 objcoord(1).y = 32.017
    objcoord(2).x = -89.204 objcoord(2).y = 32.112
    objcoord(3).x = -89.187 objcoord(3).y = 32.096

        Create Pline Into Variable myobj 0

For i = 1 to numnodes
Alter Object myobj Node Add (objcoord(i).x,objcoord(i).y)
    Next

Insert Into cables (obj) Values (myobj)

'************************************************************

MapBasic Reference says:
"The Create Pline statement creates a polyline object. If you need to create a polyline object, but do not know until run-time how many nodes the object should contain, create the object in two steps: First, use Create Pline to create an object with no nodes, and then use the Alter Object statement to add detail to the polyline object."

I'm quite new to MapBasic, so I do understand logic of this (probably some loop will be required to add node for every click), but is simply out of my current MapBasic knowledge.


texts: ArcGIS crashes with no error number in Visual Studio

texts:
I've been developing Add-ins for a while with Visual Studio 2008 and ArcGIS 10.
But in my new Add-in development process, I receive no error from Visual studio and ArcGIS crashes with error "Arcgis encountered serious application error and need to be closed".

I tried setting breakpoints and the program stops correctly at the breakpoint. But when I use F8 to proceed, I can identify the line which has the error but no error screen from Visual studio appears, as normal.

Can you help me? This way I can't debug and specify exact problem in my code


texts: Why don't new clipped or merged QGIS shapefile layers display their names in layers table?

texts:
QGIS 2.18.12 - When I use the Clip or Merge Shapefile tools in QGIS, the set of Shapefiles is created with my chosen name and in the folder I have designated.
However, in the Layers Panel, the layer is shown as 'Clipped' or 'Merged'.  Is there a good reason for this?  Wouldn't it be clearer if the Shapefile name was shown instead?  It can be confusing if I have more than one clipped or merged layer.


texts: Model script works in QGIS 2.14 but not 2.18

texts:
I created a model that worked fine in QGIS 2.14. With help I included a script in the model that prevented the model from running unless a feature was selected. But after upgrading to 2.18 the model no longer runs unless I remove the script.

Model

Error

Script

The above does exactly what I want it to do in 2.14 so I don't understand why it won't in 2.18. I have a couple of other models that use the same script for the same purpose and have the same issue. I have tried this on my PC at home and on my work PCs and still have the same issue.


texts: some fields from node and to node output are equal to 0

texts: I am generating a flow path using the flow direction grid and flow accumulation raster as input in stream to feature tool. I just want to know what it means to have 0 values in the from node and to node field and also the arcid values are 0.


texts: Map Not rendered Leaflet but renders correctly in changing the zoom level of browser

texts:
I am loading a leaflet map rendering mapbox tiles, geojson layers (initially i dont show the map) and i also have tabular links which can be seen when the page loads, the links are list of places and when I click a certain location it opens the map and zooms to the coordinates associated with it. I used following to do so

    <button id = "closeMap">Close Map</button>

The function showOnMap is executed once we click the link, the result is something like this
.
The tiles can be seen on the top left corner only.

The strange thing is when i open the browser's console or change the zoom level of browser the map comes in right place. I must have done something stupid, i am unable to find out please help me figure out.
Here is the code.


texts: SLD intervals with greater than 255 classes

texts:
I am trying to style a raster in GeoServer using a SLD file. I have 850 individual values, but only 179 color mappings (i.e., unique RGBA codes). I want to keep all 850 values, as they are meaningful, but the SLD file fails to style with more than 255 color mappings (even if < 255 distinct colors are used).

Is there a way in which I can maintain all 850 values and still style the raster with the SLD?

The work around I can see would be to create a "styling raster" that reclasses each of the 850 values to a color-mapping indicator, but this is less than desirable due to the loss of meaningful information and the redundant raster layer.

Example code to give you a better idea of what I'm talking about:

<ColorMap>
  ...
  <sld:ColorMapEntry color="#eeb241" quantity="255" label="255th value"/>
  <sld:ColorMapEntry color="#eeb241" quantity="256" label="256th value"/>
  ...
</ColorMap>

The styling from the SLD fails here. Note that the colors are the same in the above two lines (e.g., unique color number 179), but the quantities have exceeded 255.

Edit 1

The quantity to color mappings unfortunately are not organized logically, and thus I am unsure that an extended intervals type SLD will do the trick. The ColorMapEntrys look more like:

<ColorMap type="intervals">
  <ColorMapEntry color="#ffffff" quantity="3001" label="1st value"/>
  <ColorMapEntry color="#ffffff" quantity="3002" label="2nd value"/>
  <ColorMapEntry color="#0000ff" quantity="3003" label="3rd value"/>
  ... other values ...
  <ColorMapEntry color="#ffffff" quantity="3089" label="89th value"/>
  ...
</ColorMap>

So, there you can see that white is being mapped to many values, but the quantities are not organized in a logical fashion that groups the color mapping together.


texts: Android app to get altitude of places to use in field visit?

texts: I want to verify the height of the ground which is provided to me with the help of an app if available. Of course, I am aware an  app will not give me exact value but for me, it works if it shows relative height also.


texts: Point cloud - Why don't all the points show up in AGOL or ArcGIS Pro?

texts:
Whenever I upload a LAS file in ArcGIS Pro, not all the points are shown. And there are even less points when I upload it to AGOL.

Is there a maximum of points ArcGIS Pro and AGOL can show ?

If yes, what can I do to make sure ALL the points are shown?


texts: Google Earth Engine: Get average of values by month and years using Python

texts:
I am using python to download data from GEE.
I would like to get the average of the Land Surface Temperature of all the months over the years. This is what I am doing
def returnGEEAverages(satellite, subset, gemetry, name, sdate, edate):
    cover = ee.ImageCollection(satellite).filter(ee.Filter.date(str(sdate),
    str(edate))).select([subset])
    pcts = cover.reduce(ee.Reducer.mean())
    task = ee.batch.Export.image.toDrive(
                                image=pcts,
                                folder=folderName,
                                region=geometry,
                                scale=926.625433056,
                                fileNamePrefix=name,
                                crs='epsg:4326')
    task.start()
    while task.active():
      t = task.id

I am doing two loops over the years and the months
 years = np.arange(2016, 2020, 1)
 for y in years:
  for months in range(1, 13):
    days = calendar.monthrange(y, months)
    sdate = date(y, months, 1)   # start date
    edate = date(y, months, days[1])   # end date
    folderName = urban_gdf_cities['UC_NM_MN'][k]
    geometry = returnGeometry(j)
    name = str(y)+str(months)
    lstDay = returnGEEPercentiles('MODIS/006/MYD11A1',
    'LST_Day_1km', geometry, name+'LSTDay', sdate, edate)

I would like to know if there is a better way to do this. I found a similar solution but not for python.
Google earth engine:SST by month per year


texts: Can USA streets have both, name and number?

texts:
I'm working normalizing list addresses using Tiger Geocoder normalizer and a python library, usaddress E.G.

For these two I have:

5121 33RD ST N

House_no = 5121
street_no = 33

1656 HABERSHAL DR NW

House_no = 1656
street_name =  HABERSHAL

But for these one:

705 E HARRISON ST200

I'm not sure, if Harrison is name what does it means st200?


texts: Is there any library can perform polygons buffering/ inflating/ merge/ overlap at server side?

texts:
I have a set of polygon data express in lat long. I would like to perform some tasks (polygon buffering, inflating etc etc ..) in server side and output the result to client side (client is a web page... maybe use leaflet to display the data).

I am new to GIS.


texts: Getting currentfeature's color in Expression Builder of QGIS 3.xx

texts:
I want to obtain the color code of the current feature like this :
Getting currentfeature's color in Expression Builder of QGIS?
But I am in QGIS 3.14 and I dont understand how to convert the script with the new API.
I have changed
layer = QgsMapLayerRegistry.instance().mapLayersByName(layername)[0]
r = layer.rendererV2()

in
layer = QgsProject.instance().mapLayersByName(layername)[0]
r = layer.renderer()

But I have an error with
symbol = r.symbolForFeature(feature)


texts:
How can I change the css for the CartoDB tooltip?

I found this css below in the source but am having a hard time referencing it. Are there any examples of a customized, styled tooltip?

/* CartoDB tooltip control */

div.cartodb-tooltip {
    position: absolute;
    z-index: 5;
    display: none;
  }

div.cartodb-tooltip p {
    font-family: "robotoregular", Helvetica, Arial, Sans-serif;
    font-size: 15px;
    color: #333;
    text-align:left;
    text-shadow: -1px -1px 0 #FFF, 1px -1px 0 #FFF, -1px 1px 0 #FFF, 1px 1px 0 #FFF;
  }


texts: Parameters of Normal vs. Transversal cylindrical projections

texts:
I need to parameterize a cylinder that cuts the ellipsoid. One time in normal case   (eg. like Mercator) and one time in transversal case (eg. like Transversal mercator).
So, in normal case the cut cylinder is papameterized using the geographic location of the cutting circles, but in transversal case the cylinder is parameterized by setting a scalefactor. And I wonder why?


texts: I have millions of addresses that I want to geocode and I can do that in ArcGIS but it would take me endless number of hours. Is there a faster method that you can recommend?


texts: How to merge way points from many GPX files?

texts:
I have a large number of Waypoint-layers (gpx) and would like to merge them.

Is there a way to do so using QGIS?


texts:
I'm attempting to calculate land surface temperature (LST) from Landsat 7. Band 6, the thermal band for Landsat 7, comes in high gain and low gain settings. I'm calculating the temperature for Los Angeles in summer 2000 using ArcMap 10.2.2.

I know that the low gain setting (band 6L or 61) is good for scenes with high surface brightness, such as deserts. The low gain has a useful temperature range of 130-350 K. The high gain setting (band 6H or 62) is good for scenes with low surface brightness (such as those that are fully vegetated) and has a useful temperature range of 240-320 K.

I'm not sure what constitutes low and high surface brightness, relatively speaking, and whether an urban environment would count as "low" or "high" brightness. As well, I calculated LST from Landsat 5 for the same month as the one where I'm trying to calculate LST from Landsat 7, and got a high temperature above 320 Kelvin.

Due to the high temperature I got from Landsat 5 for the same scene and for close to the same day I'm leaning towards using the low gain band for Landsat 7.

Is this a correct way to think about this, or is there a better way to to choose between high-gain and low-gain for Landsat 7?


texts:
I see that other people have mentioned losing track of the layers panel, and people have rightly pointed out that it can be restored by View > Panels > Layers.

But in my install on MacOS 10.14.3, the layers panel has disappeared entirely from the list.

Same thing when I right click on the gray area above the browser. Layers panel doesn't exist.

Unless I'm missing something here – QGIS is * functionally unusable without the layers panel *. I tried

restarting qgis
trashing the preferences
restarting the machine
reinstalling an earlier version (3.4)
but no luck!


texts: Arcpy: Processing multiple files subsequently

texts:
I have an indexed dictionary containing files and want to "mosaic" (arcpy.MosaicToNewRaster_management) all files with the same index. The number of files per index range from 1 to about 30.

dic_of_files = {"index1":["file11", "file12","file13"], "index2":["file21", "file22"]}

for key in dic_of_files.iterkeys():
    arcpy.MosaicToNewRaster_management(dic_of_files[key], otherarcpyoptions)

Mosaicing 30 files in one go doesnt seem to work, I guess number is just to high. I had the idea of writing a script where python counts the number of files per index via and the loops through processing first file 1 and 2, then processing the result with file 3 and so on till it reaches the last file.

for key in dic_of_files.iterkeys():
    key_len = len(dic_of_files[key])
    if key_len < 2:
        arcpy.MosaicToNewRaster_management(dic_of_files[key], otherarcpyoptions)
    else:
        yes.. else what?

Can anybody tell me if this is a reasonable way to solve the problem? And maybe give me a hint how to script this solution?

EDIT:
ok I've gotten one step farther. I've nested a couple of loops so that python pairs up the first two items in the list, then the next two and so on (in case of an uneven item number it pairs up the last three).

The code now looks as follows (I'm omitting the arcpy function to keep it as simple as possible):

dic_of_files = {"index1": ["file10", "file11", "file12", "file13", "file14", "file15", "file16", "file17"], "index2": ["file20", "file21", "file22"]}
dic_of_pairs = {}

def addstuff(i="", v="", mydict=""):
    if i not in mydict:
        mydict[i] = []
    mydict[i].append(v)

for key, val in dic_of_files.items():
    key_len = len(val)
    if key_len < 2:
        print "do nothing"
    else:
        if key_len % 2 == 0:
            for j in range(0, (key_len - 1), 2):
                # Even number of files, process pairs all through
                output = str(str(j) + str(j + 1))
                addstuff(key, output, dic_of_pairs)
        else:
            for k in range(0, (key_len - 2), 2):
                if k < (key_len - 3):
                    # Uneven number of files. Process pairs at first..
                    output = str(str(k) + str(k + 1))
                    addstuff(key, output, dic_of_pairs)
                else:
                    # ... and process the last three files as triplets
                    output = str(str(k) + str(k + 1) + str(k + 2))
                    addstuff(key, output, dic_of_pairs)

Sadly, one round doesn't do the trick, I have to repeat this block of code on the new dictionary (here: "dic_of_pairs"). It would be very elegant if I could use this same piece of code for the second (and all succeeding) rounds, until the every index contains just one file.

Is there a way to do this? Is overwriting the original dictionary ("dic_of_files") an option, or better to pack this code into a function?


texts: QGIS 3: Set custom project variables with input dialog

texts:
I am using QGIS 3.10.0 and Windows 10. I have a project where I want the user to define custom project variables. I want to make this user friendly, I don't want the user to need to open the project settings.

So, I'm looking for a way where upon loading a project, an input dialog pops up in which the user can set custom project variables.

I've been trying to create a project macro and using PyQt5 to load an input dialog, but I keep getting error messages. I'm not a Python programmer, so I'm not even sure what I'm doing...

Could someone help me out and give me a hint on how to proceed?

EDIT
So I've added the code below to my project macro. Now whenever I launch the project, I get three seperate input dialogs. It works, but I would rather have one single dialog with three text fields. I know I have to use QDialog instead of QInputDialog, but I don't get it to work...

import sys
from PyQt5.QtWidgets import QApplication, QWidget, QInputDialog, QLineEdit
from PyQt5.QtGui import QIcon
from qgis.core import QgsExpressionContextUtils, QgsProject

project = QgsProject.instance()

class App(QWidget):

    def __init__(self):
        super().__init__()
        self.title = 'Scherm'
        self.left = 10
        self.top = 10
        self.width = 640
        self.height = 480
        self.initUI()

    def initUI(self):
        self.setWindowTitle(self.title)
        self.setGeometry(self.left, self.top, self.width, self.height)

        self.getToponiem()
        self.getOpdrachtgever()
        self.getProjectnummer()

        self.show()

    def getToponiem(self):
        text, okPressed = QInputDialog.getText(self, "Toponiem","Toponiem:", QLineEdit.Normal, "")
        if okPressed and text != '':
            QgsExpressionContextUtils.setProjectVariable(project,'Toponiem',text)

    def getOpdrachtgever(self):
        text, okPressed = QInputDialog.getText(self, "Opdrachtgever","Opdrachtgever:", QLineEdit.Normal, "")
        if okPressed and text != '':
            QgsExpressionContextUtils.setProjectVariable(project,'Opdrachtgever',text)

    def getProjectnummer(self):
        i, okPressed = QInputDialog.getInt(self, "Projectnummer","Projectnummer:", 0, 0, 1000000, 1)
        if okPressed:
            QgsExpressionContextUtils.setProjectVariable(project,'Projectnummer',i)

if __name__ == '__main__':
    app = QApplication(sys.argv)
    ex = App()
    sys.exit(app.exec_())

def openProject():
    App()

def saveProject():
    pass

def closeProject():
    pass

EDIT2
I've been trying to call an QDialog when a project opens, but it doesn't work properly. The very simple code below does get executed but the dialog window opens and immediately closes again. How can I keep the dialog open until the user presses the button? (I will later fill the dialog with more options for the user to fill in).

from PyQt5.QtWidgets import QWidget, QPushButton

def openProject():
    w = QWidget()
    w.resize(500,500)
    button = QPushButton('OK')
    button.resize(200,32)
    button.move(80, 200)
    w.show()


texts: Clipping raster by mask layer polygon produces Nodata values in QGIS

texts:
I have a raster .tif Digital Elevation Model I am trying to clip to an irregular polygon using the GDAL\Clip raster by mask layer. I ticked 'Crop the extent of target dataset to extent of cutline'.

However, the output raster has the grid cells outside the mask layer cutline converted to 'Nodata' (apart from being annoying, it also retains the large file size of the source raster). I tried setting destination nodata value to 'None', but this resulted in values of zero instead. (Please note that I do not know how to use the Commander/Python.) Any advice please?

  processing.runalg("gdalogr:cliprasterbymasklayer","Input Raster", "Clip layer.shp","",False,True,False,5,4,75,6,1,False,0,False,"","Output raster")


texts: Mmodify GPS tracker to accept two nmea inputs and gpsmarker to include a heading vector?

texts:
I use QGIS on a research vessel and wonder whether I could modify the QGIS 1.8.0 gpsinformation widget. I want to achieve two things.

Change the GPS marker to a boat shape which points in the direction of heading recieved as NMEA HDT sentence.
Add another marker from a second position feed which represents a towed underwater vehicle position. The underwater position is reported via a USBL modem which gives slant range and bearing to the target.


texts: PyCharm and ArcGIS Pro issue: Arcpy Must Run From Inside the Conda Environment,

texts:
I am getting the following error when running my python scripts in Pycharm, does anybody know how to fix it? it is very odd, because it was working just fine, then installed some updates on ArcGIS Pro and did not work anymore.

C:\Users\C5LB\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone5\python.exe C:/C5LB/Overlay_Scriptool/OverlayManyTables.py
arcpy must be run from inside the Conda environment, which ArcGIS
uses to manage the installation of Python.
You can access this environment by launching Python from the links inside the
"ArcGIS > ArcGIS Pro" Start Menu Group:
    "Python Interactive Terminal" starts a Python session.
    "Python Command Prompt" opens a Command Prompt window initialized with Conda.

Or by running the proenv.bat file in the terminal:
C:\> C:\Program Files\ArcGIS\Pro\bin\python\scripts\proenv.bat

Traceback (most recent call last):
  File "C:/C5LB/Overlay_Scriptool/OverlayManyTables.py", line 4, in <module>
    import arcpy
  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\__init__.py", line 68, in <module>
    raise ImportError("arcpy needs to run within an active ArcGIS Conda environment")
ImportError: arcpy needs to run within an active ArcGIS Conda environment


texts: Lines that touch same polygon two times in ArcGIS Desktop?

texts:
I have lines that touch neighboring polygons. But in some cases line touch same polygons two times, which I do not need.

Thus, I would like to select those line that touch one polygons two times and remove them.

Or, I would like to select lines that touch and connect different polygons.

I am using ArcGIS 10.4.1 with python, arcpy library.

I would be happy to have any suggestions to solve the case.


texts: How to delete a feature from a GeoJSON with Python

texts:
Given a GeoJSON like following:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "city": "city1"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          113.5546875,
          63.704722429433225
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "city": "city2"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          56.953125,
          41.508577297439324
        ]
      }
    }
  ]
}

How can I delete (or generate a new GeoJSON) the feature with the city property value city2 with Python?

So far I have this, though I am a little stuck in my logical right now. Would apprecaite help or suggestions.

import json
city_name = "city2"
city_geojson = [city if city['properties']['city'] != city_name else '' for city in city_geojson['features']]


texts: Customizing ArcGIS Online application template and uploading it back to ArcGIS Online?

texts:
Is it possible to customize a template for web application in ArcGIS Online?

What I mean, none of the templates 100% suits our needs so I'm wondering if I could download the code of basic viewer, let's say, add required functions and upload it back to ArcGIS Online and use it for our maps?


texts: Using Square brackets in ArcPy update cursor?

texts:
I am looking for guidance regrading the usage or meaning of the square brackets in arcpy update cursor (Win 7, ArcGIS 10.2, Python 2.7.5)
Although I have serched long and hard with various web searches, incl https://docs.python.org/2/library/index.html, I have not been able to find a descriptive answer.

So, in the example below (from post Help with Python Update Cursor syntax) I am trying to figure what the numbers within the [] actually mean?

Is it the column number?

If so, does the count start from 0 and from the left?

if(row[4] == crow.STTYPE and row[15]=="B":
    crow.FIELD1 = "HI"
    crows.updateRow(crow)
else:
    crow.FIELD1 = "BYE"
    crows.updateRow(crow)


texts:
I have a map service running on ArcGIS Server based on MXD file (ArcMap), but when I try to connect on it from QGIS I get a message:

<ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows" version="1.1.0" language="en">
<ows:Exception exceptionCode="InvalidParameterValue">
<ows:ExceptionText>Table called GEO_ATP wasn't found.</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

Any idea what is the problem?


texts: How to get a PNG image of a polygon with the corner coordinates?

texts:
How to get a PNG image of a polygon with the corner coordinates in a separate file?

I have access to ArcGIS 10.1 and QGIS 2.8.1.

I tried this:
https://gis.stackexchange.com/a/30503/50606
, but I get an uncomplete file :

0.01192239336492888

0

0

-0.01192239336492888

15.35022462480259975

-10.64971119668243738

The 2 last lines are the coordinates of the North-West corner (WGS84). I don't know what the first 4 lines represent. And I don't know how to get the 3 other coordinates.

Using the "Generate World file" in Atlas does not seem to work.


texts: Loop through geodatabases and list feature classes

texts: Trying to find a way to use Python or SQL to loop through geodatabases and print the feature classes. I'm not sure what the difference between while and for loops are; I'm brand new to scripting. A point in the right direction would be great.


texts:
I am able to create a published layers and layer groups in GeoServer based on a SQL Server store.

When I try to preview the layer in any method (e.g. OpenLayers or KML) it comes up empty. The bounds are correct (i.e. Google Earth spins to the correct location) but the geometry is not displaying in the map.

I suspect this is to do with a configuration setting for my new layer in GeoServer but I am stumped about what it might be.

Here is the layer that I am trying to get to work:

http://185.38.105.91:8080/geoserver/SCXBC/wms/kml?layers=SCXBC:uk_parliament_initial_review

Here is an example of a sample layer on the same server which does work:

http://185.38.105.91:8080/geoserver/topp/wms/kml?layers=topp:tasmania_state_boundaries

Here is a composite screenshot of the layer settings:


texts: Why is Definition Query not returning any records?

texts:
I am currently working on a mapping project and only need 8 tenements out of many to appear. I have gone into the query builder to complete this and unfortunately it will not work. Into the builder a put Tenement = 'E59/1002' AND Tenement = 'E59/1421'. It verifies but no records are found. I am not sure why this is happening as this method has always worked for me before when i need to select multiple tenements.

Any help?


texts: Error using.getInfo function within a function Google Earth Engine

texts:
I'm applying a function to map a calculation over a list of dates. In the old calculation I just manually entered each date and ran the code. As part of if that code I would use .get('year).getInfo() to obtain the year from a manually entered data.
var startDate = ee.Date('2016-10-01');
var endDate = ee.Date(startDate.advance(1,'month'));

var startyear = startDate.get('year').getInfo();
var endyear = startDate.get('year').getInfo();

I've tried do this again within my function
  var start = ee.Date(date);
  var end = ee.Date(date).advance(1,'month');
  var date_range = ee.DateRange(start,end);

  var startyear = start.get('year').getInfo();
  var endyear = start.start.get('year').getInfo();

But I'm getting the error Unknown variable references: [].
I'm guessing it's got something to do with the .getInfo function and client vs server side. However, to obtain year as an object for use, I need to use .getInfo to obtain it. Does anyone know any work arounds?
The link to my code is here


texts:
I am working in QGIS 2.18.0 version.

My analysis involves using vector filed calculator for some operation. This new filed I am creating with an output filed type as Decimal number(integer). I wrote expression and it works fine but when I save it the field type changes to Whole number.

In above image SlopeRati was Decimal type and I created one more filed SlopeRati1 and stored orginal decimal values as STRING. After saving this later doest not change but earlier filed becomes Whole number.


texts: How do I add a custom made north arrow to an ArcGIS Pro Project?

texts:
I need to be able to add a custom north arrow to achieve a standardized look for maps I make for my company.  I've been reading through the ArcGIS Pro documentation and have come to a loss on what I'm missing.  The documentation clearly states as follows.

I am attempting to add the north arrow to the "Favorites" style.  However, nowhere do I see the option for "New Item".  All I can find is the add, create, or import options.


texts:
I often have a problem when it comes to renaming featureclasses or SHPs in ArcGIS. I also have a problem when for example I have to transform 20 Feature classes or a whole dataset into a new coordinate system and want to append "_ETRS89" to the end of all the transformed feature classes.

I know there are python scripts which do this, as well variable substitution methods in model builder. But isn´t there also the option of putting this in the Tool´s dialog in ArcMap?


texts: Google Earth Engine: Error exporting a large Feature Collection to CSV

texts:
On GEE, I want to export daily weather variables for 18 years for about 30,000 rectangular polygons. That would be exporting a FeatureCollection to a csv file with ~200 million rows (30000 places * 365 days * 18 years). The export job took 2 days and then failed with this error: Error: Internal error. Is there an efficient way to export such a large feature collection?
To know if size of the export was a problem, I also tried another export with bigger rectangular polygons (which numbered about 6,000), resulting in  ~40 million rows. That job finished in about 2 hours. So I am not sure if size of the dataset is the problem here, or if something else is at play.
Here's how I am exporting:
// grids: import a FeatureCollection with 30,000 grids
// geometry: import a bounding box

var era5 = ee.ImageCollection('ECMWF/ERA5/DAILY')
                   .select(['maximum_2m_air_temperature',
                            'dewpoint_2m_temperature'])
                   .filterDate('2000-01-01', '2017-12-31')
                   .filterBounds(geometry);

var first_image = era5.first();
var scale_def = first_image.projection().nominalScale();

var era5_over_grids = era5.map(function(image) {
  return image.reduceRegions({
    collection: grids,
    reducer: ee.Reducer.mean(),
    scale: scale_def,
    crs: 'EPSG:4326'
  })});

Export.table.toDrive({
  collection: era5_over_grids,
  description: 'era5_gridded',
  fileFormat: 'CSV'
});


texts:
I have a point cloud and would like to move this to the local coordinate system for example from Revit facing toward the first quadrant. We have a file with 3 ground reference points and on-site marked with blocks. So it is possible for us to get the coordinates from the point cloud so, we have:

3 ground reference points in RDNew
3 reference points in the point cloud from the 3 blocks
an Origin point for the absolute 0,0,0 coordinates

We want to move the point cloud to the local coordinate system and rotate towards the first quadrant.

How can I perform this transformation with for example a proj4 string?


texts:
I'm a complete newb when it comes to things like Python, however I'm keen to learn more. I have some experience with ArcMap's Model Builder and have started playing around with QGIS Graphical Modeler (as we have limited licenses at my new job).

I'm trying to automate the identification of predominant land uses within subcatchments (shape file), based upon the property types i.e. predominantly detached houses, predominantly semi-detached, predominantly flats and terraced houses (this is for a hydraulic modelling application).

I've so far come up with a way to split AddressBase points based upon a house type field, into 4 separate shape files. I then plan to count the number of each type falling within each subcatchment (will result in four columns in the 'subcatchment' attribute table, with counts for each).

Next I would like to populate a fifth column, based upon the maximum value from these four counts.

So for example:

Any tips on how I could implement this?


texts: pyqgis processing, grass: --overwrite does not work

texts:
I have a QGIS plugin that tries to call the v.clean algorithm multiple times in a row and use the same output file each time the algorithm runs through. But even though I set the --overwrite parameter to true (as found in the grass manual), the output can't be overwritten and instead I get the error

 ERROR 1: A file system object called 'C:\Users\...\Desktop\output.gpkg' already exists.

in the QGIS log. I can still run the plugin, but only the snap algorithm of the following code will be applied to my line layer:

    output = self.dlg.directory.text()

    params_snap = {
        'input' : input_2,
        'type' : 1,
        'tool' : 1,
        'threshold' : 0.000095,
        'output': output,
        'error' : 'memory:',
        'GRASS_OUTPUT_TYPE_PARAMETER': 2
    }
    processing.run('grass7:v.clean', params_snap)

    params_break = {
        'input': input_2,
        'type': 1,
        'tool': 0,
        'threshold': 0.1,
        '-c': True,
        'output': output,
        'error': 'memory:',
        '--overwrite': True,
        'GRASS_OUTPUT_TYPE_PARAMETER': 2
    }
    processing.runAndLoadResults('grass7:v.clean', params_break)

The path and name of the ouput are from a QFileDialog. If the information is needed: The input file is a temporary file.

If I set a different output file for each part of the processing the code runs fine and without errors, but since I will use more than only two algorithms creating and deleting files over and over again will just slow the plugin down.

I use the current QGIS LTR 3.4.


texts:
I am trying to figure out the format of the coordinate pair below.

My task is to convert a list of these integer values into decimal degrees (below), but this is my only ground truth data point.  Does anyone have an idea of what format (it doesn't seem to be degrees minutes seconds) the above pair could be in?   The conversion should result in the decimal degrees coord pair below

?DDM?

  "current_position": {
    "latitude": 51329338,
    "longitude": -77573
  }

Decimal degrees coordinates

  "current_position": {
    "latitude": 51.548896666666664,
    "longitude": -0.12928833333333334
}


texts: QGIS plugin closing all open PyQt dialogs

texts:
In my QGIS3 plugin I need to close all open PyQt5 dialogs. The structure of my plugin is like the following code. To be clear, when I call runCloseDialogs() I want to close both dialogs opened by run() and run2() if they are open. How can I do this?

def run(self):
  self.dlg = Dialog()
  self.dlg.show()

def run2(self):
  self.dlg2 = Dialog2()
  self.dlg2.show()

def runCloseDialogs(self):
  #close all open dialogs here


texts: Is an ArcMap 10.0 Add-In compatible with ArcMap 9.3?

texts:
I'm working on an Add-In for ArcMap 10.0 in .NET version 3.5. This add-in works fine with ArcMap 10.0. Now I want to use this add-in in ArcMap 9.3. But I couldn't add Add-In .dll file in ArcMap 9.3.

QUESTION:

Is there any process to use ArcMap 10.0 Add-In in ArcMap 9.3?


texts: How to save raster results of IExtractionOp?

texts:
I'm using IExtractionOp to clip a raster by polygon:

IExtractionOp operation = (IExtractionOp)(new RasterExtractionOp());
IRasterBandCollection rasterBandCollection = (IRasterBandCollection)operation.Polygon((IGeoDataset)rasterDataset, polygon, true);

I tried saving via

// Fails with exception
IRasterDataset dataset = rasterBandCollection.Item(0).RasterDataset;
dataset.Copy("temp" + extension, (IWorkspace)rasterWorkspace);

// Successful, but does not show any pixels in ArcMap when opened
rasterBandCollection.SaveAs("temp" + extension, (IWorkspace)rasterWorkspace, rasterDataset.Format);

// Same result as SaveAs
((ITemporaryDataset)dataset).MakePermanentAs("temp" + extension, (IWorkspace)rasterWorkspace, rasterDataset.Format);

None of the three methods above worked... I can use the returned result and add it to the map, so the extraction seems to be valid, but I'm unable to store it anyhow.


texts: Large shapefile ... QGIS field calculator crashes

texts:
I have a shapefile with more than 180,000 rows. Whenever I want to apply the field calculator (e.g. for $id or $area), QGIS crashes (waiting to respond). My PC has 32GB RAM.

Any idea what I can do?


texts: How to call gdal_translate from Python code?

texts: Is it possible to use some gdal API to call gdal_translate from Python code? I do not mean simply executing the gdal_translate.exe from the file system, but rather call it somehow in code so I do not need to know the exact directory the gdal_translate executable is in?


texts: Better Distance Measurements in Web Mercator Projection

texts:
I am working with an ESRI stack, storing my layers in a sql-spatial-enabled-SDE- geodatabase (Geometry type, Web Mercator-3857).

I am building a web mapping application, so by default, the tiles are also in web mercator, 3857.

Via stored procs, i use STDistance to query distances from a user's location (Coordinates also in web mercator) to the various layers.

Problem is, due to the distortion of web mercator, my distance calcs are increasingly off, the further from the equator they are made.

I've thought of storing my layers in sql-spatial-geography (rather than geometry) type, but:

I imagine my distance queries will take much longer (distance calcs on spherical surface)
i will need to reimport a lot of data
arcgis services won't be as fast as they will need to project on the fly

If I go to Google maps, and do a distance calc, the distance returned is much more accurate, even in Nortern/southern regions, so I assume Google must be correcting for distortion caused by the web mercator projection.

My question then:  Is there a simple factor value that can be applied to distance calcs done in web mercator projection to get the "correct" distance?


texts: Dissolve Polygons by attribute leaves empty spaces

texts:
I am trying to dissolve polygons based on the same attribute using QGIS 2.8.2 (Wien).
I use Vector, Geo Processing Tools, Dissolve and then I select the attribute.
The data is Postal Sector in Uk and I need to create Postal Districts, in the attribute table I have created a field called District and I use this field to dissolve my polygons.
The problem that I have is that the new District layer has some empty spaces.
I have taken a screenshot of one of them:

The Sectors CV21 1 and CV21 2 are not part fo the new District layer and if I click on them using the "selection" icon, nothing happens.
My guess is that QGIS cannot dissolve them because CV21 2 has an intersection (exactly where its label is).
See screenshot of that area, note that those 2 lines are crossing but there is not a node:

Can anyone:
1. explain why Dissolve does not merge those polygons together
2. suggest any solution (maybe splitting CV2 1 into 2 smaller polygons and then merge them afterwards?)


texts: OpenlaRemoving FID field with WMSGetFeatureInfo function using OpenLayers?

texts:
I'm using the WMSGetFeatureInfo control of Openlayers in order to have a pupop in the map with the WMS layer info.
Automatically Openlayers inserts a FID field at the beginning of the popup.
Is there a way to hide that field?

more precisely, i'm using this function in my js map file:

   info = new OpenLayers.Control.WMSGetFeatureInfo({
        url: 'http://localhost:8080/geoserver/wms',
        title: 'Identify features by clicking',
        queryVisible: true,
        eventListeners: {
            getfeatureinfo: function(event) {
                map.addPopup(new OpenLayers.Popup.FramedCloud(
                    "chicken",
                    map.getLonLatFromPixel(event.xy),
                    null,
                    event.text,
                    null,
                    true
                ));
            }
        }
    });
    map.addControl(info);

Do you have any suggestion to remove the fid field?


texts: JMapPane doesn't update after adding new layer

texts:
After adding a new layer it doesn't show inside JMapPane, is there a way to force this?
I already tried using MapLayerEvent

styleSld = Stil.createFromSLD(inputStream);
layer = new FeatureLayer(pointCollection, styleSld);
layer.setTitle("G_POINT");

getMapPane().getMapContent().addLayer(layer);
MapLayerEvent mple = new MapLayerEvent(layer, MapLayerEvent.DATA_CHANGED);
MapLayerListEvent mplle = new MapLayerListEvent(getMapPane().getMapContent(), layer,
    getMapPane().getMapContent().layers().indexOf(layer), mple);
getMapPane().layerChanged(mplle);

The same code is used for refreshing content when drawing new feature and it works normaly.

I also tried running new ResetAction and nothing happens. If I interacted with a map by (zoom change or by clicking ResetButton or Drag) it will render the layer.

EDIT

The problem is caused by,

getMapPane().getMapContent().getViewport().setBounds(new ReferencedEnvelope(-100, 100, -100, 100,
                getMapPane().getMapContent().getCoordinateReferenceSystem()));

after I remove this it works nicely, but then need to go through the code and manually create bounds if a return is null, but it will be much easier with this part in the code.


texts: Write raster with its attribute table and symbology

texts:
I am classifying lots of rasters and I'd like to automate the legend creation and symbology displaying for easy manipulation.

Here is the kind of raster I'm dealing with:

library(raster)
r <- raster(nrow=10, ncol=10)
r[] = 1
r[51:100] = 3
r[3:6, 1:5] = 5
r <- ratify(r)

rat <- levels(r)[[1]]
rat$landcover <- c('Pine', 'Oak', 'Meadow')
rat$class <- c('A1', 'B2', 'C3')
levels(r) <- rat
r

levelplot(r, col.regions=c('palegreen', 'midnightblue', 'indianred1'))

How could I embed this symbology in the raster when writing it?


texts: Sencha / Extjs 4 with GeoExt 2 - How to use

texts:
I'm working on Sencha Architect 2 - extjs 4
i found out that GeoExt doesn't work with extjs 4. but i really need to use it so i searched and end up here https://github.com/geoext/geoext2

I'm trying to something like examples so i just add :

Ext.Loader.setConfig({
    enabled: true,
    disableCaching: false,
    paths: {
        GeoExt: "../geoext-2/src/GeoExt",
        // for dev use
        Ext: "http://cdn.sencha.io/ext-4.1.0-gpl/src"
        // for build purpose
        //Ext: "extjs-4.1.0/src"
    }
});

to my app.js file. but i keep getting "GeoExt is not defined" error message.

anyone could help me with that ?


texts:
I would like to display information of a line feature in parallel. For example, a river line feature with its different qualities. I want to avoid a parallel copy of line features. Is it possible to achieve this with display symbology options in ArcMap?
For clearance:

from: http://wbv-untere-warnow-kueste.de/fileadmin/processed/csm_6_5_GWzustand_karte_76134b7b82.jpg


texts:
I just started worrking with OpenLayers and hava a problem with markers that don´t appear.
I have no idea where is the problem.

<script type="text/javascript">

function init() {
    var center = new OpenLayers.LonLat(17.55, 53.7).transform('EPSG:4326', 'EPSG:900913');
    var map = new OpenLayers.Map("map");
    var osm = new OpenLayers.Layer.OSM();
    map.addLayer(osm);
    map.setCenter(center, 12);

    var markers = new OpenLayers.Layer.Markers("Markers");
    map.addLayer(markers);
    var size = new OpenLayers.Size(21, 25);
    var offset = new OpenLayers.Pixel(-(size.w / 2), -size.h);
    var icon = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png', size, offset);
    var table;

    console.log(map);

$.getJSON('url',

function (data) {
    console.log(data);
    $.each(data.response.venues, function (i, venues) {

        markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(venues.location.lng, venues.location.lat).transform(new OpenLayers.Projection("EPSG:4326"),new OpenLayers.Projection("EPSG:900913")), icon));
        markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(venues.location.lng, venues.location.lat).transform(new OpenLayers.Projection("EPSG:4326"),new OpenLayers.Projection("EPSG:900913")), icon.clone()));
    });
});
}

1 UPDATE: map.addLayers(markers) changed to map.addLayer(markers)
2 UPDATE: put the getJSON call in the init() function


texts: How to change a widget type of a field by using python

texts:
I'm working on a plugin for QGIS 3.x and I'm looking for an example how I can change the widget type by using Python code. The change should be applied on fields of a vector layer.

At the end it should look like that:

I found a solution which works:

setup=QgsEditorWidgetSetup( 'ValueMap', {  'map': {u'Description 1':
u'value1',  u'Description 2': u'value2'} } )
layer.setEditorWidgetSetup(3, setup)


texts: Stopping new points being added by Create Features of ArcGIS Pro?

texts:
I am using Create Features from the Edit tab of ArcGIS Pro 2.1.2, and am using the Point construction tool to add multiple points (not multipoints).

When I want to stop adding points I expected to be able to click F2 or right-click and choose Finish.

The help on Create point and multipoint features is a little ambiguous on this but when I want to stop adding points I expected to be able to click F2 or right-click and choose Finish.

However, when I click F2 nothing happens (I still seem to be in the middle of adding new points) and right-clicking on or away from any of the points just added does not offer Finish in the one item context menu.

I think this behaviour may have changed from the last ArcGIS Pro version that I used to add points but, in any event, is there an obvious way to stop creating new point features that I am overlooking?


texts: OSM roadways shapefile - Completing missing road names

texts:
I have a shapefile for the road network in Portugal that I downloaded from a OpenStreetMap mirror site (geofabrik), dated about a month ago.
A little context beforehand for those not familiar with the road naming system in Portugal -- since it may not be the same for your country--: We have municipal roads, which usually belong within city limits; national roads, that serve a purpose of connecting the different cities, without using highways; primary and complementary itineraries, which act as fast roads --like highways, but with a few differences--; and highways. There are others, but these are the main ones.
Currently I'm using the "type" field (motorway, primary, secondary,...), to establish some sort of hierarchy, but this doesn't work for me because I need to separate itineraries from highways, and national from municipal.
What I would like to know is if there is any way for me to complete the "name" (XX Street, Y Avenue,...) or the "ref" (Highway 1, Municipal Street 273,...) fields, since they are lacking a lot of values, and thus being able to use the Query Builder to sort them all out.
I searched a lot before writing this question and most of the posts I see are about missing roads, and this is not the case: as far as I know, my road network is complete.
I am also aware that there is a new shapefile available to download on the geofabrik, dated yesterday, but unless they had a huge overhaul, the road names will still be missing.


texts:
I am not able to understand how to extract values of NTL and Population for computing NLDI.
I have the NTL and Population grid both in same projection and resolution.
How should I proceed in QGIS?


texts: Aerial Photograph Footprint Size Calculation

texts:
I was wondering what an easy formula to calculate the dimensions of the ground covered by an aerial photograph.

I know the following:

camera sensor size (mm x mm)
camera focal length (mm)
height above ground (m)

How can I combine these variables to calculate metres by metres covered by each individual image captured?


texts:
I've got a PostGIS database that's been loaded with a couple of the continent downloads from geofabrik.  One of these continents is Antartica - the pbf download is 29 MB but the size in the database is around 1.5 GB (50 times larger) - at this rate the whole world would be over 2 TB.  Is this normal?  Is there any way to reduce the size?

One thing got me wondering - when I look at the nodes table, there are over 4.5 million nodes, which includes not just the, approx 1600, node features but also every vertex of every way in the continent too.  Is this normal?

All the tags are being extracted into an hstore.

The process was done using osmosis.


texts: Geotagging photos with smartphone camera?

texts:
I have a motorola turbo 2 smartphone. I am exploring options for automatically geotagging photos for site visits. The goal is to add these photos into Google Earth such that I can see where they were taken. Here is the information I have collected so far:

1) The format that I need these photos to be in is EXIF format

2) A Google Photos feature exists, but it doesn't seem to include options for geotagging (you can only geotag one photo at a time, for social media sharing purposes)

3)Other apps exist (Photomap, OSMTracker, GPS Map Camera, Photo Exif Editor) but none of them seem to provide a way to do this.

My question is: Am I on the right track to figure out a solution? Should I direct my efforts elsewhere?


texts: Editing file geodatabase without using tools from Esri?

texts:
I am looking for a tool (Python library etc.) that would allow me to edit a file geodatabase without Esri software?

It is possible?


texts: Is it possible to use the ArcObjects Runtime SDK on a Windows Store Application (windows 8, 8.1)

texts: I'm looking to develop a simple windows 8 tablet application with arcobjects. Is it possible to use the ArcObjects Runtime SDK? If not, is there an ArcObjects SDK for Windows store Applications?


texts: Build web app engine from google earth engine

texts:
I want to make a web app engine using code from google earth engine which will find the landsat8 image based on filter date on HTML.

I use a python program to deploy the web app.

Any ideas how to make it?


texts: Is it possible to map a leaflet provider (precipitation) on top of a polygon in an interactive map?

texts:
I've just started using the leaflet package in R to create interactive maps and I'm wondering if it's possible to plot one of the leaflet providers on top of a plotted polygon rather than under it as a basemap. Here's some example code:

library(raster)
library(leaflet)

ph <- getData("GADM", country = "PHL", level = 1)

leaflet(data = ph) %>%
  addProviderTiles("OpenWeatherMap.Precipitation") %>%
  addPolygons(fillColor = "red",
              fillOpacity = 1,
              color = "#FFFFFF",
              weight = 1.3)

This simply plots a shapefile of the Philippines (in my real example, it's a choropleth map). The precipitation basemap was found here (toward the bottom of the list; OpenWeatherMap.Precipitation). Here's what the output looks like:

It might be a little difficult to see, but the precip layer is plotted as a basemap, so all of that info is going under the ph shapefile. I'm trying to find a way so that the precipitation layer goes OVER the Philippines shapefile. I'm not sure if this is even possible, but it seems like I should be able to have a plain OSM basemap with the shapefile plotted on top and the precipitation layer on top of the .shp. As it is now, I know that I can alter the opacity of the shapefile, but this isn't much of a solution for a detailed choropleth map. Is this possible, and if so, how do I do it?


texts: Automatic four-colour (political boundary) rendering with PostGIS

texts: This is probably a long shot, but is there a way to automatically select colours for neighbouring regions using the four colour theorem, in PostGIS? The data I'm using is suburb boundaries in OpenStreetMap (boundary=administrative, admin_level=10).


texts: qgis2web can't create map with WMS in EPSG:3857

texts:
I have a simple vector  data (EPSG:2177) and I want to add a WMS server:

http://mapy.geoportal.gov.pl/wss/service/img/guest/ORTO/MapServer/WMSServer

I choose a EPSG:3857 CRS, and still I have a blank page when updating my map in qgis2web plugin (sometimes it works on Leaflet but i need it on OpenLayers 3).

Any clues what am I doing wrong?

Update:

When I open WMS with Qgis I see 8 available WMS (also EPSG:3857):


texts:
I'm new to Mapinfo and Mapbasic as well.
I need to build a tool like in ArcMap (pythonaddins).
The goal is to have a gui application which has three  dropdown lists.
First is a list of layers in workspace.
Second one is a set of fields inside the selected table
This
Third is a list of values that have to be written in a specific field (according to selected features).

I already discovered that Mapbasic can make dialog windows but things that I tested were working for one time while I needed a tool to be like an instrument which will stay always in Mapinfo window.
Any suggestions?


texts: In leaflet does the addControl automatically handle events via the leaflet API?

texts:
I need to add additional overlay layers to be controlled by the controls within the leaflet map boundaries, but the code utilizes a leaflet-custom.js file. I am assuming it handles the event and triggers creating DOM elements for the overlay but I cannot verify as the custom file is a minified file so I don't know what configurations were made.

When I add my wms tile layers to the controls the wms tile layers do not display as a result of 403 errors with the Forbidden message, the urls are all modified, and the type changes from images/pngs/jpegs to text/html.

Therefore is there any way to programmatically add geolocation with with the events?


texts: Dissolving complex polygons using PyQGIS?

texts:
I have a workflow to eliminate small polygons within a layer and merge them with adjacent polygons. My input is the result of a polygonization of a raster (has no overlaps) and it has one field "class".

Input:

In the first step I loop over all features and if a polygon is of class 0 (blue), and is too small, i change its class to a suitable neighbour.

Next I dissolve the features, according to their changed class values.

As you can see, the blue areas are dissolved successfully, however the island in the middle appears to have changed value. Actually there are now two polygons laying on top of each other there, a light green one, and the original one. The light green basically filled a hole during the dissolve process. This happens only for a couple of polygons, and not all of them.

Here's the code to change the values of the polygons:

def KFB_for_1_class (input, field_, c, area_threshold):
    feats = []
    dis_temp = r"z:\dis.shp"

    layer = iface.addVectorLayer(input,"","ogr")
    layer.startEditing()

    layer_cpy = iface.addVectorLayer(input,"","ogr")

    print "cleaning class "+str(c)
    for f in layer.getFeatures():
        if f[field_] == c and f["area"] < area_threshold:
            org_int = f[field_]
            candidates = []
            for f2 in layer_cpy.getFeatures():
                if f.geometry().touches(f2.geometry()):
                    if f2[field_] not in candidates and f2[field_] != org_int:
                        candidates.append(f2[field_])

            candidates.sort()
            if len(candidates) >0:
                f[field_] = candidates[0]
            for i in (candidates):
                if abs(org_int - i) <= abs(f[field_] - org_int):
                    f[field_] = i
                    f["area"] = area_threshold + 1

            layer.updateFeature(f)
            feats.append(f["id"])

            print "feature " +str(f["id"]) + " changed from " + str(org_int) + " to " + str(f["int"])

    print feats

    layer.commitChanges()

To dissolve I use this:

processing.runalg("qgis:dissolve", input, False , "class", dis_temp)

What I tried:

using other dissolve tools within QGis (OGR, SAGA, [GRASS I could not get to work])
cleaning Topology with v.clean according to Dissolve polygons in QGIS (GEOS) with tolerance

Both did not work. The different tools all ended up creating the same error.

How can I resolve this?


texts: PostGIS layer visible in QGIS Browser panel but displayed as not connected in DBManager

texts:
I have added my database as a workshop and you can see the layers in public schema. but whereas in Db manager it is saying NOT CONNECTED. how can I resolve this?
this is running on localhost. on windows with all latest version software


texts: problem display tile when config tilecache for WMS

texts:
i have a problem with config tilecache for WMS
when i have shown it. but there are proplem when display tiles.

any idea?

this is my config for tilecache.cfg

[world]

type=WMS
url=myserver/mapserv.exe?mode=map&map=C:\Mapffiles\countries.main.map
layers=countries
bbox=-180, -90, 180, 90
srs=EPSG:4326
extension=png
reproject=true
levels=15
size=256,256
metaSize=5,5
sendfile=yes
metaTile=true

when i display it in QGIS then it fails as above


texts: Moving geometry at same distance in different coordinate systems

texts:
I haven't been making tools in ArcPy for a while and stuck in one case.

My goal is to create a serie of geometries in a center of current extent by running a script. I.e when I move ArcMap's extent and press a custom tool button, four points at certain distances between each other and one polyline are added to specific layers in a center of extent. This is made for saving time when cartographer has to digitize same shapes multiple times.

The algorithm seems to be pretty simple:

Find the center coordinate of current extent
Create multiple points changing only X coordinate of each new point by adding float numbers of distances in meters
Create line between first and last points
Add geometries to layers.

As a result a plain horizontal line with four points should appear in layers and ArcMap's workspace. Something like that:

But other users can set another coordinate system in dataframe and in this case a center of extent can be a coordinate containing decimal degrees or meters. Script can add fixed decimal values to an degree-X coordinate of extent's center and finally distances can be completely different while I need a strict metrical distances in any coordinate system.

I tried to transform dataframe's extent with .projectAs() method and got a center coordinate in Equidistant projection. Then performed calculations for center point, generated geometries and projected them to inital dataframe's coordinate system with a same .projectAs() method.

Unfortunately if I change dataframe's settings by pointing another coordinate systems and run script the result will be different every time. I think there should be an easier solution rather then transforming geometries.

Here is a piece of code:

import arcpy

def get_center_coordinate(dataframe):
    ex = dataframe.extent.projectAs(arcpy.SpatialReference(check_crs))
    x_center = (ex.XMax + ex.XMin) / 2
    y_center = (ex.YMax + ex.YMin) / 2

    return (x_center, y_center)

# current dataframe data
mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
init_sr = df.spatialReference

# metric crs
check_crs = 54027
check_crs_obj = arcpy.SpatialReference(check_crs)

# calculation of new coordinates and creating geometry objects
center_cor = get_center_coordinate(df)
pnt_leftmost = center_cor[0] - 1.145
pnt_gas = arcpy.PointGeometry(arcpy.Point(pnt_leftmost, center_cor[1]), check_crs_obj).projectAs(init_sr)
pnt_arm_first = arcpy.PointGeometry(arcpy.Point(pnt_leftmost + 0.29, center_cor[1]), check_crs_obj).projectAs(init_sr)
pnt_arm_second = arcpy.PointGeometry(arcpy.Point(pnt_leftmost + 0.59, center_cor[1]), check_crs_obj).projectAs(init_sr)
pnt_cross = arcpy.PointGeometry(arcpy.Point(pnt_leftmost + 2.29, center_cor[1]), check_crs_obj).projectAs(init_sr)
list_line_new = arcpy.Array([arcpy.Point(pnt_leftmost, center_cor[1]), arcpy.Point(pnt_leftmost + 2.29, center_cor[1])])
polyline_segment = arcpy.Polyline(list_line_new, check_crs_obj).projectAs(init_sr)

# making insert cursors
ic_gas = arcpy.da.InsertCursor('gas_pnt','SHAPE@')
ic_arm = arcpy.da.InsertCursor('arm_pnt','SHAPE@')
ic_cross = arcpy.da.InsertCursor('cross_pnt','SHAPE@')
ic_segments = arcpy.da.InsertCursor('segments','SHAPE@')

# adding geomtries to layers
ic_gas.insertRow([pnt_gas])
ic_arm.insertRow([pnt_arm_first])
ic_arm.insertRow([pnt_arm_second])
ic_cross.insertRow([pnt_cross])
ic_segments.insertRow([polyline_segment])

del ic_gas, ic_arm, ic_cross, ic_segments


texts:
Im creating a plugin for VBS2. In the program I can output a grid in either LAT LONG,UTM or MGRS. I need to be able to convert to BNG. I have managed to create a TKinter application in python that works using Proj.4 but now need to create it in C++ as a DLL.

The LatLong im using (51.20650N 1.81906W) is a known point, the BNG conversion is approx SU 127 452.

#include <proj_api.h>

double x = 51.20650; //atof(GetStrArgument(input, 0));
double y = 1.81906;  //atof(GetStrArgument(input, 1));

char *pj_latlongc = "+proj=longlat +datum=WGS84";
char *pj_UTMc = "+proj=utm +zone=29 +ellps=WGS84";
char *osc = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000    +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs";

projPJ pj_OS = pj_init_plus(osc);
projPJ pj_UTM = pj_init_plus(pj_UTMc);
projPJ pj_latlong = pj_init_plus(pj_latlongc);

x *= DEG_TO_RAD;
y *= DEG_TO_RAD;

int p = pj_transform(pj_latlong, pj_OS, 1, 1, &x, &y, NULL);

Unfortunately the result is completely wrong and doesn't make sense.Could someone help on this?


texts: How to get Netherlands neighborhood names based on Addresses?

texts: CBS provides key figures on the neighborhood level for the whole country(NL). I want to incorporate this data for analysis but my datapoints include addresses and postcodes and latitudes&longitudes, not neighborhood names. Is there an easy way to match each address to the corresponding neighborhood?


texts: Simplify raster size to plot using ggplot geom_raster()

texts:
I want to create a map using 'ggplot2'. My data contains shapefile (points data), polygons (European countries). Underneath the vector data, I wanted to place the raster showing extends of forests over Europe, available here. The raster size of disk is 150 MB, has dimension of 46000 and 65000. To use it in ggplot, I need to first convert it to dataframe: raster_df<-as.data.frame(raster).
But, I got an error: Error: cannot allocate vector of size 11.1 G. I do not understand why I got this error? How to make sure that I can plot my raster with vector data using ggplot2?
Dummy example:
# read raster
forest <- raster(paste(getwd(), "eea_r_3035_100_m_forest-area-2015_p_2015_v1_r1.tif", sep = "/"))

# convert to dataframe for ggplot2()
forest_df <- as.data.frame(forest, xy = TRUE)

# Plot the raster data on top of European countries
ggplot(data = europe) +
geom_fs() +
geom_raster(data = forest)

EDIT:
I have increased my memory.limit(12000) and this allows me to create a dataframe or a stars object that I can both plot with ggplot2. But, the display still takes forewer. When I visualize the map in ArcGIS or QGIS, the plotting takes a second. Isn't there something as efficient to plot rasters and use them in a map in R? tmap::tm_raster() seems promissing but other suggestions are welcome!


texts: Creating a new layer from a spreadsheet and an existing shapefile in QGIS

texts: I have a map of US zip codes and a spreadsheet that associates data with specific zip codes. I want to show the zip codes where the most data exists using a gradient, but I am not even sure how to get it onto the map in the first place.


texts:
Running QGIS (using currently 3.14.1, but this happens in older versions as well) on Windows 10, I sometimes encounter the problem of a task (algorithm) running endlessly without the possibility to stop it. The only workaround is to force QGIS to close, losing temporary layers.
In this specific case, I started vectorize - it's running for two hours and I don't know for how long it will continue. I would like to change parameters and start it again. However, the Cancel button is not active.

There are other cases of tasks that I canceled, but they still continue to run in the background. So it seems to be a general issue of QGIS.
Is there a solution to cancel running algorithms without closing the program?


texts: Extract the values from a raster based on polygons and add those values directly to the polygon shapefile table in ArcGIS?

texts: I understand that the ZonalStatistics (ArcToolbox -> Spatial Analyst Tools->Zonal->Zonal Statistics) in ArcGIS (10.3) allows the extraction of raster values within zones defined by a polygon shapefile (takes the average over the specified areas).  ZonalStatistics creates a new raster file based on the polygons.  I'd like to directly add the raster values to my polygon shapefile as a new column.


texts:
We are developing a tracking software and we are the stage where we test the reliability of our location stigma on the following web-maps

Google Map
Nokia HERE
Open Street Map (OSM)

We observe that our location appears accurate on Google maps, however when switched to Nokia or Openstreet there is a discrepancy of around 40-50 metres from our actual location.
Most importantly our GPS position accuracy is around 2m.

Could it be inaccuracy of map calibration or any similar problems known by anyone?


texts: Is the order in which ArcMap loads Extensions, Toolbars and Commands defined?

texts: When you extend Esri's ArcMap with both an Extension, a Toolbar and a number of Commands, is the order in which these categories are loaded by ArcMap well-defined? (I.e. are Commands always loaded before Toolbars, Extension always loaded before Commands, etc.)


texts: Points dataset legend to polygon dataset

texts:
I'm working with ArcMap 10.5.1
I have a point dataset in a gdb, with a particular legend. The color in the legend represent the area of effect that intersects that point. I also have a polygon dataset, in the same gdb, that represents said areas of effect. The coincident field with the area of effect ID is present in both datasets.
In my layout I'd like to showcase the points with a different color depending on the area of effect they belong to, and then, below this layer, the areas of effect dataset with the same colors but some transparency. The colors were randomly generated in the symbology options.
I exported my point dataset legend to a lyr file to transfer it to the polygons, but it does not work because of the different types of geometry. And entering the colors manually would be too much work.
Is there any way I can transfer the legend from one dataset to the other?


texts: Showing only uppercase text as labels in QGIS?

texts:
Using QGIS 2.14 and see there is now rule based labelling. I wish to show labels of towns that have capital letters (eg CRAWLEY) and not show the labels of towns with lowercase letters (eg Faygate)

The attribute table does not distinguish the different towns apart from their name. Is there a way of doing this ?


texts:
The answer to this question tells us how to wrap labels on spaces. Is it possible to wrap labels on defined spaces? For example, wrap on just the first space or on just the second space.

Secondly, is there a way of controlling the amount of vertical space inserted when wrapping labels? Wrapping on a character (or on a control-character) inserts too much vertical space.

Nick.


texts: Showing attribute table in form using PyQGIS

texts:
Is there any way to show the attribute table of a layer in a QTextBrowser or a QTableView on form?

I did this but, (shows attribute table in a separate dialog) I have no idea, how to show it in my form.

iface.showAttributeTable(iface.activeLayer)


texts: Why do simplify polygon results differ between ArcGIS for Desktop and PostGIS?

texts:
I am mainly interested in why and how the Simplify Polygons tool in ArcGIS for Desktop handles the simplification process, including preventing and resolving topological errors better than the equivalent ST_SimplifyPreserveTopology in PostGIS destroys topology when trying to simplify to similar levels.

Use case scenario: PostGIS database queried for vectorr GeoJSON polygon overlay in Leaflet. One can pass the ST_SimplifyPreserveTopology to the geometry field but my tests have shown that preprocessing data in ArcGIS gives me smaller datasets and topological correctness.

Here is an example:

Case 1 - Original data

Case 2 - ST_SimplifyPreserveTopology(the_geom,0.001)the_geom

note: querying from CartoDB

Further, when simplifying Case 1's data in ArcGIS for Desktop, there are no topology errors produced (ArcGIS for Desktop can resolve all errors).


texts: Polar projection polygon visualization issue

texts:
I am using Natural Earth coastline polygon (NE_10m_land), loaded from a shapefile in QGIS. All displays well (in, e.g., unprojected WGS84, EPSG 4326) - see below:

But when visualising in a polar projection (here, EPSG:3571, WGS84 in a North Pole LAEA Bering Sea projection), QGIS seems to lose the ability to distinguish between the inside and outside of polygons.

This does not seem to be a problem with the shapefile ( not well formed polygons) as I can visualize it using ArcGIS just fine.

So this seems to be  a problem with the rendering engine in QGIS. Any hints?

Edit:
It appears that it is indeed caused by Antarctica. Manual removal of Souther Hemisphere helps.


texts: Calculating Slope of Street Segments using ArcGIS Desktop?

texts:
Using a DEM, I need to assign slope to individual segments in a street centerline feature class. It will be used to account for increased travel time for steep slopes in a network analysis.  I am using ArcGIS 10.0, ArcInfo license, with 3D Analyst and Spatial Analyst extensions available.

My initial attempt was to use the "Add Surface Information" tool.  However, the resulting average slope values seem suspect in some areas, for example if a stream crosses the street segment.  I would like to try assigning elevations to street segment endpoints then calculating slope.  I have read about this method Calculating average slope along line using ArcGIS Desktop?) but need further details on the process.


texts:
I need to calculate the PIPEID field to equal the G_KEY_COMP field if PIPEID is null or PIPEID and G_KEY_COMP do not match.  I think the first condition of the PIPEID being null is covered by the second condition of the fields not matching?  I'm a novice on this topic and can't get the syntax I have to work. Here's what I've piecemealed together:

Pre-Logic Script Code:

def update(PIPEID,G_KEY_COMP):
 if PIPEID != G_KEY_COMP
  return G_KEY_COMP
else:
 return PIPEID

PIPEID = update(!PIPEID!, !G_KEY_COMP!)

I'm posting a screenshot of the field calculator window since I'm not able to get the indentation right in the codeblock.


texts: Spatial join between polygon with multipolygon in QGIS

texts:
I have two layers in QGIS - wgs84 representing the study area and nycIncome which I would like to clip based on where polygons intersect those of wgs84

Both layers are in the same CRS and have almost identical geometry, but they are not the same.  Whenever I try to run a spatial join using FIRST in the MMQGIS plugin, however, QGIS crashes.

Is there a reason why this is happening and what would be the correct way to take the intersection between a polygon and multipolygon?

.


texts: PC and Mac treating raster brick input differently in R

texts:
All - this is likely an idiotic question, but, I've been working on developing some GIS code in R that interacts with an online source.  So, sprinkled in various places, I have

myBrick <- brick("http://somesite.com/image.jpg", crs=proj)

This works well and good. And then...I moved to a PC. And chaos ensued.  The issue is that brick() on a PC appends additional text to the filename. So, I get

Error in .local(.Object, ...):
`C:\path\to\my\working\director\http://somesite.com/image.jpg` does not exist in the file system

This is.... odd. Same identical code. Entering strings. Never referencing working directory.

What is going on here? What is brick() doing differently? How can I stop the madness?


texts: Setting interactivity on sublayer created in CartoDB?

texts:
I am trying to set the interactivity on a sublayer that I created on CartoDB web interface so that my events and tooltips receive the correct fields as arguments of the callback functions. Here is how I have tried to do this:

    layers[1].setInteraction(true)
    poverty_layer = layers[1].getSubLayer(0)
    poverty_layer.set("interactivity", ["cartodb_id, namelsad10, hh_median"])

    poverty_layer.on "featureOver", (e,ll,pos,data)->
      console.log data

I assumed that by doing this I would see three data fields –– cartodb_id, namelsad10, and  hh_median –– but I only see cartodb_id. I already managed to set interactivity successfully on sublayers that are created on the clientside. This is the one use case where it hasn't been successful.


texts:
Please can someone show me how to export a file from R in the format of a shapefile in order to be used in GeosegregationAnalyzer (which requires a shapefile).

I have merged attribute data with the shapefile on R as shown below and want to export this new file called "Classmap01"

Brismap01 = st_read("lsoa01/lsoa01.shp")

Class01 = read.csv("ClassData2001.csv")

Classmaps01 = merge(Brismap01, Class01, by.x="lsoa01cd", by.y="code")


texts: How to troubleshoot "RuntimeError: The field is not nullable [Column_Name]"?

texts:
I have a python (arcpy) script that extracts data from a database table, processes it, and inserts into a feature class. It works in my development environment but not in production. I recently removed a "IS NOT NULL" constraint from the database table so now the column is allowed to accept NULL values.

I verified the constraint was removed from this table. The column is nullable. And then I checked one more time... just to be sure

Where do I look next? How do I troubleshoot this issue?

Update 1 - NULL values are okay. I'm not trying to prevent NULL from hitting the database. If a value doesn't exist for a field, I want NULL.

Update 2 - The field in the database is set to allow NULL values. I verified it again. When viewing the FC properties in the catalog, within ArcMap, the field "Allow NULL values" is set to No. I'm guessing that's the problem. Why are the two out of sync and what do I do to fix it?


texts: Creating a Line feature in PyQGIS QGIS 3.x

texts:
In QGIS GUI, a line feature is created from this icon

1.jpg

How to Create a line with a fixed distance from a vector layer in python console.

This my example:

from qgis.core import *

from qgis.PyQt.QtCore import QVariant

capa = QgsVectorLayer("Point?crs=epsg:32718", "temp", "memory")

pr = capa.dataProvider()
pr.addAttributes([
                  QgsField("codigo",  QVariant.Int),
                  QgsField("nombre", QVariant.String)])
capa.updateFields()

..

In Toolbars / mDigitizeToolBar the is action mActionAddFeatures this create a line feature, i need how to call from Python Console, with parameters .


texts: Replacing more than one escape character when labeling in ArcMap using VBscript?

texts:
I have written a label expression to match font color to the symbology for a point file based on unique values. Problem is, occasionally the comments will include one or the other special characters ("&" or "<"). Sometimes there can be hundreds of points per map that may or may not contain those special characters.

I have replaced one of the escape characters using the snippet:

function findLabel([FORESTER],[NAME]):
NewString = Replace( [NAME] ,"<","&lt;")
if([FORESTER] = "AG") then
    findLabel = "<CLR red= '211' green='255' blue='190'>"&NewString&"</CLR>"
elseif ([FORESTER] = "RS") then
    findLabel = "<CLR red= '76' green='230' blue='0'>"&NewString&"</CLR>"

But how do I write the code to replace 2 values in the same variable?


texts: ArcGIS display multiple popup at the same time?

texts:
I'm using ArcGIS API v4 and the document said ArcGIS API allows view can only have one instance of popup at a time

All Views contain a default popup. This popup can display generic content, which is set in its title and content properties. When content is set directly on the Popup instance it is not tied to a specific feature or layer.

I've been searching for similar questions but still got no clearly answers? Is there a way to implement opening multiple popup on ArcGIS map at the same time?


texts: Stereo Compilation Software with PostGIS

texts:
Does anyone know of stereo compilation software that will allow for the writing into a PostGIS database?

Both Datum and Stereo Analyst will write to an ESRI Geodatabase.  Standard stereo compilation software will write to a DGN, but I am unaware of anything that will write directly into a PostGIS database.


texts:
I have researching possibility to fetch raster file from WMTS server via GDAL. I have found some old post about that GDAL doesn't support WMTS (WMS only). However, I also find this information: http://gdal.org/frmt_wmts.html, where is mentioned that

  [...] Since GDAL 2.1, access to WMTS layers is possible with the GDAL WMTS client driver (needs Curl support). [...]

However I couldn't find any source of this GDAL 2.1, newest version on gdal.org is 2.0.

Can you tell my where can I find GDAL 2.1 or how to enable WMTS support for older version of GDAL? Now i use GDAL 1.11.2 on OS X.


texts: How do you create a query that will result with the adjacent parcels?

texts:
I have a shapefile with parcel data, and I am trying to create a query that will render a list or something similar that will list adjacent parcels to a specific parcel.

I want to be able to select the my parcel and the query to result with all the adjacent property owners.

I am using ArcGIS 10.1


texts: Is it necessary to do an accuracy assessment for an unsupervised classification done for a Landsat 5 image?


texts:
Can you tell me how to create a circle in leaflet which can be dragged across the map?

I do not want to use Circlemarkers as i want to specify radius in meters and also the size  of circle should increase as we zoom in the map according to the distance it covers.

I want the radius in meters to show how many markers lie within that circle or that region. i want the dragging functionality so that user can reposition the circle across the map.

So basically i need to create a circle on the map which can be dragged across the map

i am using leaflet with mapbox for creating this application.


texts: Get colortable from three-band RGB raster (brick) with R

texts:
I have a three band RGB Raster Geotiff. I will provide the sample data below, but I'm guessing my problem should be solvable generically. I would like to turn this three band raster into a single band raster with a color table.

The original three bands can each range from 1 to 255, resulting in over 16M (256^3) theoretically possible combinations of Red, Green and Blue. Of these 16M+ Combinations, little more than 11'000 combinations are actually used. Be that as it may, the resulting color table allows only 256 values in total.

Is there an algorithm to reduce the 16M+ possible RGB-values to an appropriate set of 256 values?

I'm guessing that I'd first have to decide on a colour palette with 256 colours and then somehow assign each RGB Value to one of these colours. Before I reinvent the wheel, I thought I'd ask here whether this problem already has been solved. I'd like to use R, the package raster and any other necessary package for this.

Note: Interestingly, the producers of said map also provide a single band geotiff with just 13 (!!) different colours in the colortable, but it doesn't look much different from the three band Raster (see below).

This is the three band RGB raster:

This is the single band raster with a color table:

See the sample data on the swisstopo homepage. In case the Link breaks: https://shop.swisstopo.admin.ch > "Digital national maps > Swiss Map Raster > Swiss Raster Map 50 > Sample Data > Swiss Map Raster 50, Sample data. The three band RGB Raster is in the folder SMR_50KREL_254dpi_uncompressed/SMR50_LV95_KREL_10L_Mosaic.tif. The single band Raster with a colortable is here: SMR_50KOMB_508dpi_LZW/SMR50_LV95_KOMB_Mosaic.tif


texts:
I am involved in a bigger university mobile application project and my task is to design and develop a mobile module regarding map and location-based services in a native Android application.

Shortly my part of app should display different kind of info on a georeferenced map (Google Map or Openstreetmap) based on Oracle data:

Shortest path to reach a certain event from a certain location
Different layers with domain information to the event (like parking places around it, bus lines ecc. )
GetFeatureInfo option

I am new in this kind of topic and I would like to ask you an architectural question regarding it.

I thought about using Geoserver to define layers, styles, datasets etc., and WMS protocol to generate the map over a website using  HTML 5, OpenLayers, and jQuery-mobile.

But how can I accomplish this in a native android app?

Building a webpage where I display the map using HTML 5, OpenLayers and jQuery-mobile to interact with Geoserver and then using the Webviews in Android to display the webpage in my app? Would be this a possible solution? Possible drawbacks?

What are the alternatives using Geoserver?


texts: May I change rgrass7 location in R script

texts:
I'm using Grass-GIS 7.2 in R with rgrass7. I would like to use rasters from 2 different location. Is there a way to do that without brutality ?
At this time it works with :

library(rgrass7)
# initialization and the use of my first location
initGRASS(gisBase = "/usr/lib/grass72", home = tempdir(),
          gisDbase = "/home/user/Documents/futurSahel/grassDB/",
          location = "landsat8", mapset = "PERMANENT")

# list of all raster loaded in GRASS
execGRASS("g.list", parameters = list(type = "raster"))

ndvi.landsat <- readRAST("ndvi.rescal",  close_OK=TRUE, drivername="GTiff")

##Using the 2nd location
initGRASS(gisBase = "/usr/lib/grass72", home = tempdir(),
          gisDbase = "/home/user/Documents/futurSahel/grassDB/",
          location = "spot", mapset = "delaye", override=TRUE)

# list of all raster loaded in GRASS
execGRASS("g.list", parameters = list(type = "raster"))

ndvi.landsat <- readRAST("SPOT6_P_20151113",  close_OK=TRUE, drivername="GTiff")


texts:
I have a shapefile with polylines representing rivers all over a region.
In the attached image, there is an example of one river.

The longest route is marked by the two red dots on this river, but it is divided into many polylines especially along this route. I would like to find the total length of this river (and many many others) without having to manually delete all tributaries then merge the longest route.

Is there an efficient method for accomplishing this task? Is it possible to do this task to many rivers at the same time?

I'm using ArcGIS 10.1 with an ArcInfo license. I have access to VBA, but not Visual Studio.


texts: How do I create a polygon from a buffered point?

texts:
I have a point and buffered it using st_buffer. How do I create a polygon geometry from that buffered point? The result of st_buffer is the polygon I need to save for future use.

I will be placing the new geometry in the geometry column of an existing table (created specifically for this purpose) and use the id of the point as the foreign key to link it back to the original point. Potentially, there may be multiple buffers for each point so point's id will not be the rows PK.

The destination table will look something like this...

CREATE TABLE geom_buffered_pts
(
  oid bigserial NOT NULL,
  point_id bigserial NOT NULL,
  geom geometry,
  CONSTRAINT geom_buffered_pts_pk PRIMARY KEY (oid)
)


texts: How to connect to a Portal for ArcGIS js API

texts:
I have some issues connecting to my ArcGis Portal JavaScript API. I don't know what to specify in:

require(["in here"],
    function (in here) {

});
And also where to put authentification parameters ?
they make it look quite simple in that Esri documentation: https://developers.arcgis.com/javascript/3/jsapi/portal-amd.html.

And I tried many times without success.
Can anybody help me with that ?


texts: spTransform error (unable to find inherited method for function spTransform)

texts:
I am trying to reproject a shapefile.  For example, the land shapefile from Natural Earth.

library(sf)
shape <- read_sf(".../ne_10m_land.shp")

shape_proj<-spTransform(shape, CRS("+proj=gnom +lat_0=90 +lon_0=-50"))

I keep getting the following error, even when I try with CRSobj= or using other ways to express the CRS.

  Error in (function (classes, fdef, mtable)  :    unable to find an
  inherited method for function ‘spTransform’ for signature ‘"sf",
  "CRS"’


texts: Dynamically update points with xy coordinates in table

texts:
I have an Oracle database table with columns for 'X' and 'Y' coordinates. Now I want to display these points in ArcMap. I can use the 'Display XY data' tool for that.

My table in the database with XY coordinates will be updated every day, so I have to update my points in ArcMap with the new coordinates, too.

Is there any way to automatically update my points in ArcMap based on my table in database?


texts: Re-project GeoJSON from pixel-space to lat-long in Python

texts: I have created a GeoJSON of linestrings in pixel space using Python.  However, I now need to re-project the file to match the same projections and extent of the .tiff file I wish to layer the GeoJSON on top of.  Is there a straightforward way to do this using rasterio, geopandas, or GDAL?  I would rather not have to loop over the coordinates of every linestring since that is time consuming and computationally expensive.  How should I go about this?


texts: How to correctly append data from ArcInfo GENERATE file to PostGIS table using ogr2ogr

texts:
Using CentOS 8.0 with PostgreSQL 12, PostGIS 3, and GDAL/OGR 3.0.2. I am trying to append data in an ArcInfo GENERATE file to an exisiting table.

My ArcInfo GENERATE file input.gen looks like this:

1 11.250000 58.282526
2 -168.750000 58.282526
3 -78.750000 69.094843
4 -123.750000 35.264390
5 -78.750000 31.717474
6 -33.750000 35.264390
7 -57.844843 -0.000000
8 -20.467474 -0.000000
9 11.250000 20.905157
10 11.250000 -20.905157
11 42.967474 -0.000000
12 56.250000 35.264390
13 80.344843 -0.000000
14 101.250000 31.717474
15 101.250000 69.094843
16 146.250000 35.264390
17 -137.032526 -0.000000
18 -99.655157 0.000000
19 -123.750000 -35.264390
20 -78.750000 -31.717474
21 -33.750000 -35.264390
22 -78.750000 -69.094843
23 11.250000 -58.282526
24 56.250000 -35.264390
25 101.250000 -69.094843
26 101.250000 -31.717474
27 122.155157 0.000000
28 146.250000 -35.264390
29 159.532526 -0.000000
30 -168.750000 20.905157
31 -168.750000 -20.905157
32 -168.750000 -58.282526
END

The PostGIS table is created with:

CREATE TABLE myschema.points
(
    gid bigint NOT NULL
        CONSTRAINT points_pkey
            PRIMARY KEY,
    geog geography not null
);

CREATE INDEX points_gist
    ON myschema.points (geog);

Constraints:

The table column names cannot be changed.
The ID in the ArcInfo GENERATE file must be retained and loaded into the gid
If the gid already exists, ignore the error and continue. Similar to PG's ON CONFLICT DO NOTHING

The command I tried:

ogr2ogr -doo "ACTIVE_SCHEMA=myschema" \
        -append \
        --config PGSQL_OGR_FID gid \
        -nln points \
        -sql "SELECT ID AS gid FROM input" \
        -f "PostgreSQL" PG:"dbname= user= password= host= port=5432" \
        input.gen

But this results in a not-null constraint violation

  ERROR 1: ERROR:  null value in column "gid" violates not-null constraint

  DETAIL:  Failing row contains (null, 0101000020E610000000000000008026408522DDCF29244D40).

  ERROR 1: INSERT command for new feature failed.

  ERROR:  null value in column "gid" violates not-null constraint

  DETAIL:  Failing row contains (null, 0101000020E610000000000000008026408522DDCF29244D40).

  Command: INSERT INTO "points" ("geog" ) VALUES ('0101000020E610000000000000008026408522DDCF29244D40'::GEOGRAPHY)

  ERROR 1: Unable to write feature 0 from layer centers.

  ERROR 1: Terminating translation prematurely after failed
  translation from sql statement.

Note:

It seems the ID is not properly selected, so I tried id, fid, ogc_fid, and oid in the -sql statement instead, but no luck.
Then I tried -preserve_fid, because it sounds like this option will preserve the fid of the input data. That command actually goes through without an error, but it does not load the data as in the input, instead it starts at 0.

How do I correctly select and insert the ID in the inputfile?

The ogr2ogr docs and ArcInfo GENERATE have helped me.

This may sound like my previous question, but it is not. It is a different issue.


texts:
I was hoping that the following command would give me an alpha-channel PNG which was transparent where the DEM has "nodata":

gdaldem color-relief -alpha -of png dem.tif palette.dat out.png

I get an RGBA PNG alright, but the alpha channel is completely opaque! Why?

I am reasonably convinced that there are "nodata" pixels in dem.tif, because if I add an nv entry to the palette (with some garish colour), I get a lot of this colour in my output, and right where I expect it to be.

gdaldem --version says GDAL 2.2.1, released 2017/06/23.


texts: Problem While Selecting from table using ST_Transform